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Abstract 


The  contents  of  this  report  covers:  (i)  development  of 
optimal  geometry  for  crowned  helical  gears;  (ii)  method  for  their 
generation;  (iii)  tooth  contact  analysis  (TCA)  computer  programs 
for  the  analysis  of  meshing  and  bearing  contact  of  the  crowned 
helical  gears  and  (iv)  modelling  and  simulation  of  gear  shaft 
deflection. 

The  developed  method  for  synthesis  is  used  for  determination 
of  optimal  geometry  for  crowned  helical  pinion  surface  and  is 
directed  to  localize  the  bearing  contact  and  guarantee  the 
favorable  shape  and  low  level  of  the  transmission  errors. 

Two  new  methods  for  generation  of  the  crowned  helical  pinion 
surface  have  been  proposed.  One  is  based  on  application  of  the 
tool  with  a  surface  of  revolution  that  slightly  deviates  from  a 
regular  cone  surface.  The  tool  can  be  used  as  a  grinding  wheel 
or  as  a  shaver.  Other  is  based  on  crowning  pinion  tooth  surface 
with  predesigned  transmission  errors.  The  pinion  tooth  surface 
can  be  generated  by  a  computer  controlled  automatic  grinding 
machine . 

The  TCA  program  simulates  the  meshing  and  bearing  contact  of 
the  misaligned  gears.  The  transmission  errors  are  also 
determined . 

The  gear  shaft  deformation  has  been  modelled  and 
investigated.  It  has  been  found  the  deflection  of  gear  shafts 
has  the  same  effects  as  those  of  gear  misalignment. 
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SUMMARY 


The  tope'  gy  of  several  types  of  modified  surfaces  for 
helical  gers  are  proposed.  The  modified  surfaces  allow 
absorption  of  linear  or  almost  linear  function  of  transmission 
errors  caused  by  gear  misalignment  and  deflection  of  shaft. 
These  surfaces  result  in  the  improved  contact  of  gear  tooth 
surfaces.  The  principles  and  corresponding  programs  for  computer 
aided  simulation  of  meshing  and  contact  of  gears  have  been 
developed.  The  results  of  this  investigation  are  illustrated 
with  numerical  examples. 

1.  INTRODUCTION 

Traditional  methods  for  generation  of  involute  helical  gears 
with  parallel  axes  provide  developed  ruled  tooth  surfaces  for  the 
gear  teeth  (Fig.  1.1).  The  tooth  surfaces  contact  each  other  at 
every  instant  along  a  line,  L,  that  is  the  tangent  to  the  helix 
on  the  base  cylinder.  The  surface  normals  along  L  do  not  change 
their  orientation.  The  disadvantage  of  regular  helical  gears  is 
that  they  are  very  sensitive  to  misalignments  such  as  the 
crossing  or  intersection  of  gear  axes.  The  misaligned  gears 
transform  rotation  with  a  linear  function  of  transmission  errors 
(a  main  source  of  noise)  and  the  bearing  contact  is  shifted  to 
the  edge  of  the  teeth.  The  frequency  of  transmission  errors 
coincides  with  the  frequency  of  ‘-ooth  meshing.  The  actual 
contact  ratio  (the  average  number  of  teeth  being  in  mesh  at  every 
instant)  is  close  to  one  and  is  far  from  the  expected  value. 
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These  are  the  reasons  why  we  have  to  reconsider  the 
canonical  ideas  on  involute  helical  gears  and  modify  their  tooth 
surfaces.  Crowning  of  the  gear  surfaces  is  needed  to  negate  the 
effects  of  transmission  errors  and  the  shift  of  contact  between 
the  gear  tooth  surfaces.  Deviations  of  screw  involute  gear  tooth 
surfaces  to  provide  a  new  topology  that  can  reduce  the  gear 
sensitivity  to  misalignment  will  be  developed.  Theoretically, 
the  modified  tooth  surfaces  will  be  in  contact  at  every  instant 
at  a  point  instead  of  a  line.  Actually,  due  to  the  load  applied 
between  meshing  teeth,  the  contact  will  be  spread  over  an 
elliptical  area  whose  dimensions  may  be  controlled.  Methods  for 
gear  tooth  surface  generation  that  provide  the  desirable  surface 
deviation  are  proposed.  For  economical  reasons  only  the  pinion 
tooth  surface  is  modified  while  the  gear  surface  is  kept  as  a 
regular  screw  involute  surface. 

2.  BASIC  CONCEPTS  AND  CONSIDERATIONS 
2.1  Simulation  of  Meshing 

The  investigation  of  influence  of  gear  misalignment  requires 
a  numerical  solution  for  the  simulation  of  meshing  and  contact  of 
gear  tooth  surfaces.  The  basic  ideas  of  this  method  (Litvin, 
1968)  are  as  follows: 

(1)  The  meshing  of  gear  tooth  surfaces  is  considered  in  a 
fixed  coordinate  system,  Sf.  Usually,  the  generated  gear  tooth 
surfaces  may  be  represented  in  a  three  parametric  form  with  an 
additional  relation  between  these  parameters  -  Gaussian 
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coordinates.  Such  a  parametric  form  is  the  result  of 
representation  of  a  gear  tooth  surface  as  an  envelope  of  the 
family  of  the  tool  surface  (the  generating  surface)  and  two  from 
the  three  Gaussian  coordinates  are  inherited  from  the  tool 
surface . 

The  continuous  f  angency  of  gear  tooth  surfaces  is 
represented  by  the  following  equations 


.(1) 


2) 


r  (u^,0^,i|»^,<t>^)  =  £  (2.1.1) 


‘  2  '  v  2 '  *  2 '  *  2  ‘ 


n  (Ujf0^fip^r<})^)  —  n  ^  ‘^(u 2/02' ^2' ^2^' 


n(1>  =  n(2) 


(2.1.2) 


f 6^U 1' 9 1' ^ 1^  =  °  (2.1.3) 

f ? (u 2*  0  2 »  ^2^  =  ®  (2.1.4) 


Here:  u^  and  0^  are  the  tool  surface  curvilinear 
coordinates,  iJk  is  the  parameter  of  motion  in  the  process  of 
generation  of  the  gear  tooth  surface,  is  the  angle  of  rotation 
of  the  gear  being  in  mesh  with  the  mating  gear. 

Equations  (2.1.1)  to  (2.1.4)  provide  that  the  position 
vectors  r ^ ^  and  r^2^  and  surface  unit  normals  n ^ ^  and  n^2^  ace 

iV  l-W  A, 

equal  for  the  gear  tooth  surfaces  in  contact  (Fig.  2.1).  Vector 
equations  (2.1.1)  and  (2.1.2)  yield  five  independent  equations 
and  the  total  equation  system  is 
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fi(Ui,0i,<j>i,U2»02,<f)2,^l'^2^  =  i  •  ®  *  [  1  f  5  ]  / 


f e ( u  1 ' e  i  r  ^  i )  =  0,  f 7(u2» 9 2'^2^  ”  °  (2.1.5) 

An  instantaneous  point  contact  instead  of  a  line  contact  is 
guaranteed  if  the  Jacobian  differs  from  zero,  i.e.  if 

D(f,  -  rf. if. -  fEj-tf.-)) 

— i 1—5 *  o  (2.1.6) 

D(u^,0  ]/'f’i'u2'92'^2,'*,2 

If  the  inequality  equation  (2.1.6)  is  observed,  then  the  system 
of  equations  (2.1.5)  may  be  solved  in  the  neighborhood  of  the 
contact  point  by  functions 

ul(<t>i>  ,u2(4>i)»^i<^i)  ^  •  »  •  ♦  .»4>2^i^ 

These  functions  are  of  class  (at  least  they  have  continuous 

derivatives  of  the  first  order).  Functions  (2.1.7)  and  equations 
(2.1.5)  enable  calculation  of  the  transmission  errors  (deviation 
of  2  ( ^  1  ^  ^rom  t!ie  Prescribed  linear  function)  and  the  path  of 
the  contact  point  over  the  gear  tooth  surface. 

For  the  case  when  the  gear  tooth  surface  is  a  regular  screw 
involute  surface,  it  may  be  directly  represented  in  a  two- 
parametric  form  and  the  number  of  equations  in  system  (2.1.5)  may 
be  reduced  to  s4  . 
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2.2  Simulation  of  Contact 


IXie  to  the  elastic  approach  of  the  gear  tooth  surfaces  their 
contact  is  spread  over  an  elliptical  area.  It  is  assumed  that 
the  magnitude  of  the  elastic  approach  is  known  from  experiments 
or  may  be  predicted.  Knowing  in  addition  the  principal 
curvatures  and  directions  for  two  contacting  surfaces  at  their 
point  of  contact  we  may  determine  the  dimensions  and  orientation 
of  the  contact  ellipse  (Litvin,  1968). 

The  determination  of  principal  curvature  and  directions  for 
a  surface  represented  in  a  three-parametric  form  is  a  complicated 
computational  problem.  A  substantial  simplification  of  this 
problem  may  be  achieved  using  the  relations  between  principal 
curvatures  and  directions,  and  the  parameters  of  motion  for  two 
surfaces  being  in  contact  at  a  line.  One  of  the  contacting 
surfaces  is  the  tool  surface  and  the  other  is  the  generated 
surface . 


Helical  gears  with  modified  gear  tooth  surfaces  will  be 
designed  as  surfaces  being  in  point  contact  at  every  instant. 
The  point  of  contact  traces  out  on  the  surface  a  spatial  curve 
(the  path  of  contact)  whose  location  must  be  controlled.  The 


tangent  to  the  path  of  contact  and  the  derivative  of  the  gear 


ratio  —  (m».  (<J>,))  may  be  controlled  by  using  the  relationship 

a<p  zi  i 

between  principal  curvatures  and  directions  for  two  surfaces  that 


are  in  point  contact  (ref.  2).  Here: 


^2 

“1 


fux) 
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is  the  gear  ratio. 


and  they  are  also  represented  by  a  piece-wise  linear  function 
(Fig.  2.2(b)).  Transmission  errors  of  this  type  cause  a 
discontinuity  of  the  gear  angular  velocity  at  transfer  points  and 
vibration  becomes  inevitable.  The  new  topology  of  gear  tooth 
surfaces  proposed  in  this  report  allows  the  absorption  of  a 
linear  function  of  transmission  errors  that  results  in  a  reduced 
level  of  vibration.  This  is  based  on  the  possibility  to  absorb  a 
linear  function  by  a  parabolic  function. 

Consider  the  interaction  of  a  parabolic  function  given  by 


(2.3.3) 
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with  a  linear  function  represented  by 


The  resulting  function 

2 

A<)>2  =  -  a<j>  ^  (2.3.5) 


may  be  represented  in  a  new  coordinate  system  (Fig.  2.3): 


*2  = 


(2.3.6) 


where 


+  2  =  A*2 


(2.3.7) 


/  i  v  2 

VJe  consider  that  A<t>2  '  =  —a 4>  is  a  predesigned  function  that 
exists  even  if  misalignments  do  not  appear.  The  absorption  of 
function  A  4>  2  ^  ^  =  t>4>  by  the  parabolic  function  means 

that  gear  misalignment  does  not  change  the  predesigned  parabolic 
function  of  transmission  errors.  Thus  the  resulting  function  of 
transmission  errors  a  4>  2  =  ^2^  +  A  ^ 2 ^  wiii  keep  its  shape  as  a 
parabolic  function  although  the  gears  are  misaligned.  The 
resulting  function  of  transmission  errors  may  be  obtained 

by  translation  of  the  parabolic  function  a^^. 

The  absorption  of  a  linear  function  of  transmission  errors 
by  a  parabolic  function  is  accompanied  by  the  change  of  transfer 
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points.  The  transfer  points  determine  the  positions  of  the  gears 

where  one  pair  of  teeth  is  rotating  out  of  mesh  and  the  next  pair 

is  coming  into  mesh.  The  change  of  transfer  points  is  determined 
b  b 2 

with  A  *  ^  =  and  A 4> ^  =  — t*ie  cycle  of  meshing  of  one  pair 

2ir 

of  teeth  is  4>  ^  =  —  i  e  1,2.  It  may  happen  that  the  absorption 

of  a  linear  function  by  a  parabolic  function  is  accompanied  with 

a  change  that  is  too  large.  If  this  occurs  the  transfer  points 
and  the  resulting  parabolic  function  of  transmission 

errors,  i|;  (tji  ),  will  be  represented  as  a  discontinuous  function 

for  one  cycle  of  meshing  (Fig.  2.4).  To  avoid  this,  it  is 
necessary  to  limit  the  tolerances  for  gear  misalignment. 


2.4  Misalignment  of  Regular  Helical  Gears 

Regular  helical  pinion  and  gear  can  be  represented  by  their 
surface  position  vectors  and  normal  vectors  in  coordinate  system 


and  S2  as: 


2 

COS  * 


cosU  -  *  )  [t  sins  +  r.*  ]  +  r.sin< 


r  COS*  ~~~'Tp  vc,  ^p~..pp 

X1  2C, 

COS  * 

tri]  =  Yl=  FoT*^  sinUP  '  *c)[tPsinBP  +  ri*P]  +  ricos*P 

Z1  c  .  2a 
l  2  sin  6  2 

tp(cos8p  *  sin  *n  >  +  rl*psi"  *nVo 

1  p 

(2.4.1) 

[cos*  cosB  cos*  +  sin*  sin*  1 
.  .  n  P  P  n  vp 


[n, ]  =  n.  !  =  -cos*  cosB  sin*  +  sin*  cos 


n  p  Tp 


rP  TP 


+  cos*  smB 
I  n  p 


(2.4.2) 


8 


[r2] 


2 

^  2 


2, 

COS  dj 

_ _r 

COS  ip 


sin(4>G  -  4>c )  [-tGsinB  +  r24>G]  +  r2s:''n< 


n 


2 

cos  ^ 

cos*  sinUG  "  *cX-tGsin8p  +  r 2^]  +  r2cos*G 

C  •  2 

2  sin  Bd 

t^(cos8  +  sinz*  - ^—2.)  -  r 

G  p  yn  cosg 


■24.Gsin%ntgep 


(2.4.3) 


[n2]  = 


2x 

n0 

2y 

Ln2z 


cos*  cos8  cos*.  +  sin*  sin*,, 
n  p  G  yn  G 


-cos*  cose  sin*„  +  sin*  cos*_ 
Yn  p  G  yn  G 


-cos*  sine 
n  p 


(2.4.4) 


where  *^  and  *^  are  the  pressure  angles  in  gear  tooth  normal 

section  and  transverse  section  respectively;  Bp  is  the  helical 

angle  at  the  pitch  cylinder  of  the  pinion  and  the  gear;  rj  and  r2 

are  the  radii  pitch  cylinder  of  the  pinion  and  the  gear 

respectively;  *^,  and  t^  are  the  surface  parameters  of  the 

pinion  tooth  surface;  *  and  t„  are  the  surface  parameters  of  the 

G  G 

year  tooth  surface. 

When  the  helical  pinion  and  gear  are  in  mesh,  with  their 
axes  misaligned,  their  position  vectors  and  normal  vectors  can  be 
transformed  to  the  fixed  coordinate  system  S^.  The  basic  ideas 
have  already  been  discussed  in  the  Section  2.1.  And  the  real 
approach  of  transformation  and  the  matrices  to  describe  the 
misalignment  and  gear  rotation  will  be  given  in  Section  3.8. 

It  is  found  that  when  the  misalignment  occurs  the  regular 
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pinion  and  gear  surface  cannot  contact  in  tangency.  That  is,  in 
Sf,  their  normals  can  not  be  equal  in  all  circumstances.  In  this 
case,  only  the  gear  tooth  edge  contacts  pinion  tooth  surface  in 
tangency.  Therefore,  in  the  fixed  coordinate  system,  there  are 
four  equations  to  describe  the  contact,  that  is,  the  equalities 
of  three  position  vectors  describing  gear  tooth  edge  and  pinion 
tooth  surface  as  well  as  the  zero  product  of  gear  tooth  edge 
tangency  and  pinion  tooth  normal. 

The  computer  aided  simulation  of  meshing  of  misaligned 
helical  gears  with  regular  tooth  surfaces  shows  that  the 
transformation  of  rotation  is  accompanied  with  large  transmission 
errors.  There  are  two  sub-cycles  of  meshing  during  the  complete 
meshing  cycle  for  one  pair  of  teeth.  These  sub-cycles  correspond 
to  the  meshing  of  (1)  a  curve  with  a  surface,  and  (2)  a  point 
with  the  surface.  The  curve  is  the  involute  curve  at  the  edge  of 
the  tooth  of  the  gear  and  the  point  is  the  tip  of  the  gear  tooth 
edge.  The  transmission  errors  for  the  period  of  a  cycle  are 
represented  by  two  linear  functions  (Fig.  2.5).  The 
transformation  of  rotation  will  be  accompanied  with  a  jump  of  the 
angular  velocity  of  the  driven  gear  and  therefore  vibrations  are 
inevitable . 

The  results  of  computation  are  presented  for  the  following 
case.  Given:  number  of  teeth  N-^  =  20,  N2  =  40,  diametral  pitch 
in  normal  section  Pn  =  10  in-*,  gear  tooth  length  L  =  10/Pn  the 
helical  angle  6^  =  15°,  the  normal  pressure  angle  ^  =  20°.  The 
gear  axes  are  crossed  and  form  an  angle  ay  =  5  arc  minutes.  The 
computed  transmission  errors  are  represented  in  table  2.4.1. 
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TABLE  2.4.1  -  TRANSMISSION  ERRORS  OP  REGULAR  HELICAL 
GEARS  WITH  CROSSED  AXES 


41  T  • 
deg 

-5 

- j - 

- 3— 

7 

Aij>^, 

sec 

4.90 

3.06 

1.22 

-0.61 

-2.45 

-4.29 

2.5  Surface  Deviation  by  the  Change  of  Pinion  Lead 


Helical  gears  in  this  case  are  designed  as  helical  gears 
with  crossed  axes.  The  crossing  angle  is  chosen  with  respect  to 
the  expected  tolerances  of  the  gear  misalignment  (Ay  is  in  the 
range  of  10  to  15  arc  minutes).  The  gear  ratio  for  helical  gears 
with  crossed  axes  may  be  represented  (Litvin,  1968)  as 


M 


rb2  Sin  Xb2 
rbl  Sin  Xbl 


(2.5.1) 


where  rbi  and  xbi  are  the  radius  of  the  base  cylinder  and  the 
lead  angle  on  this  cylinder,  i  e  1,2.  j ~  Xpl|  =  Ay* 
Here:  X  ^  is  the  lead  angle  on  the  pitch  cylinder.  The  advantage 
of  application  of  crossed  helical  gears  is  that  the  gear  ratio  is 
not  changed  by  the  misalignment  (by  the  change  of  Ay).  The  tooth 
surfaces  contact  each  other  at  a  point  during  meshing.  The 
disadvantage  of  this  type  of  surface  deviation  is  that  location 
of  the  bearing  contact  of  the  gears  is  very  sensitive  to  gear 
misalignment.  A  slight  change  of  the  crossing  angle  causes 
shifting  of  the  contact  to  the  edge  of  the  tooth  (Fig.  2.6). 

The  discussed  type  of  surface  deviation  is  reasonable  to 
apply  for  manufacturing  of  expensive  reducers  of  large  dimensions 
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«•-  A _ « 


when  the  lead  of  the  pinion  can  be  adjusted  by  regrinding.  While 
changing  by  regrinding  the  parameters  r^  and  t^ie 
requirement  that  the  product  r^  sinA^  must  be  kept  constant. 
Then,  the  gear  ratio  M21  will  be  of  the  prescribed  value  and 
transmission  errors  caused  by  the  crossing  of  axes  will  be  zero. 

Theoretically,  transmission  errors  are  inevitable  if  the 
axes  of  crossed  helical  gears  become  intersected.  Actually,  if 
gear  misalignment  is  of  the  range  of  5  to  10  arc  minutes,  the 
transmission  errors  are  very  small  and  may  be  neglected.  The 
main  problem  for  this  type  of  misalignment  is  again  the  shift  of 
the  bearing  contact  to  the  edge  (Fig.  2.6). 


3.  GENERATION  OF  PINION  TOOTH  SURFACE  BY  A  SURFACE  OF  REVOLUTION 
3.1  Basic  Consideration 

The  purpose  of  this  method  for  deviation  of  the  pinion  tooth 
surface  is  to  reduce  the  sensitivity  of  the  gears  to 
misalignment.  Also  the  transmission  error  must  be  kept  to  a  low 
level  and  stabilize  the  bearing  contact.  This  investigation 
shows  that  this  goal  may  be  achieved  by  the  proposed  method  of 
crowning  but  the  bearing  contact  cannot  cover  the  whole 
surface.  The  reason  for  this  is  that  the  instantaneous  contact 
ellipse  moves  across  but  not  along  the  surface  (Fig.  3.1). 

The  proposed  method  for  generation  is  based  on  the  following 
considerations.  It  is  well  known  that  the  generation  of  a 
helical  gear  may  be  performed  by  an  imaginary  rack-cutter  with 
skew  teeth  whose  normal  section  represents  a  regular  rack-cutter 
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for  spur  gears  (Fig.  3.2(a)).  We  may  imagine  that  two  generating 

surfaces,  EQ  and  E^,  are  applied  to  generate  the  gear  tooth 

surface  and  the  pinion  tooth  surface,  respectively  (Fig. 

3.2(b)).  Surface  E  is  a  plane  (a  regular  rack-cutter  surface), 

and  E  is  a  cone  surface.  Surfaces  E_  and  r  are  rigidly 

p  G  p 

connected  and  perform  translational  motion,  while  the  pinion  and 
the  gear  rotate  about  their  axes  (Fig.  3.3).  The  generated 
pinion  and  gear  will  be  in  point  contact  and  transform  rotation 
with  the  prescribed  linear  function  ^  2 ( ^ 1 ) *  HoW6ver'  due  to  9ear 
misalignment,  function  <p 2 ( 4>  1 )  becomes  a  piecewise  linear  function 
(Fig.  2.2(a))  that  is  not  acceptable.  To  absorb  a  linear 
function  of  transmission  errors  (Fig.  2.2(b)),  a  predesigned 
parabolic  function  of  transmission  errors  is  used.  For  this 
reason  a  surface  of  revolution  that  slightly  deviates  from  the 
cone  surface  is  proposed  (Fig.  3.2(c)).  The  radius  of  the 
surface  of  revolution  in  its  axial  section  determines  the  level 
of  the  predesigned  parabolic  function.  The  pinion  crowning 
process  may  be  accomplished  by  grinding,  shaving  or  lapping. 

3.2  Principle  of  Generation  and  Used  Coordinate  Systems 

Consider  two  rigid  connecting  surfaces  E„  and  E  .  The 

G  p 

generating  surface  E0  is  a  plane  and  generates  the  helical  gear 

tooth  surface  that  is  an  involute  screw  surface.  Surface  E  is  a 

P 

surface  of  revolution.  Initially,  we  consider  that  E  is  a  cone 

P 

surface  and  E^  and  eq  contact  each  other  along  a  straight  line 
that  is  the  generatrix  of  the  cone. 
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Fig.  3.4  shows  the  generating  surfaces  E„  and  E  .  Fig.  3.5 

G  P 

illustrates  the  process  for  generation.  While  the  rigidly 
connected  generating  surfaces  perform  a  translational  motion,  the 
pinion  and  gear  rotate  about  their  axes  0-^  and  0 2 , 
respectively.  The  parameter  of  motion  of  the  cutter,  S,  and  the 
angles  of  rotation  of  the  pinion  and  the  gear,  <f>^  and  are 

related  as  follows: 


(3.2.1) 


where  r^  and  r2  are  the  great  centrodes  radii,  the  cutter 

centrode  is  the  straight  line  that  is  tangent  to  the  gear 

centrodes.  Point  I  is  the  instantaneous  center  of  rotation. 

Coordinate  systems  and  S2  are  rigidly  connected  to  the  helical 

pinion  and  the  helical  gear,  whereas  coordinate  systems  Sc  and  Sf 

are  rigidly  connected  to  the  tool  surface  and  fixed  frame.  The 

generating  surface  1  (a  plane)  and  E  are  covered  with  a  set  of 

G  p 

contact  lines  LG2  and  Lp^  respectively the  instantaneous  lines 

of  tangency  of  surfaces  EG  and  e2  and  E  and  E.  (Fig.  2.3, 
a,b).  The  location  of  these  lines  depends  on  the  value  of 


parametric  <j>  and  <k  (<k  and  $  are  related).  Line  is  the 

P  P  G  op 

line  of  tangency  of  generating  surfaces  E~  and  E  . 

G  p 

When  e  ^  and  E2  are  generated,  at  any  moment,  one  point  on  line 
LpG  are  generating  corresponding  point  Mpi  and  MGi  on  the  helical 
gear  surface  E  ^  and  E2*  When  the  £  ^  and  i'2  are  meshing  without 
misalignment  and  MGi  contact  each  other  in  turn.  Now  it  is 

clear  why  LpG  is  not  parallel  with  the  edge  of  Eq.  The  reason  is 
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that  in  this  way,  the  contact  ratio  of  crowned  helical  pinion  and 
regular  helical  gear  will  be  higher. 


3.3  Tool  Surface 


The  pinion  generating  surface  is  a  cone  and  may  be 
represented  in  an  auxiliary  coordinate  system  (Fig.  3.7)  as 

follows : 


r 


d 


f  u  cosesina 
? 


d-u  cosa 
P 


-u  sin9  sina 
P  P 


1 


(3.3.1) 


where  0  <  u  <  -  ,  0  <  e  <  2n .  The  surface  normal  is 

cosa 

represented  by 


N 


d 


u  sina 
P 


fcosacose 

sina 

cosasine 


P 

P. 


The  unit  surface  normal  is  (provided  u^sina  *  0) 


(3.3.2) 


'cosacose 

P 

n.  =  sina 
~d 

cosasine 

L  PJ 


(3.3.3) 


Figure  3.8  illustrates  the  installment  of  the  conic  tool  in 
coordinate  system  Sc  step  by  step  as  follows: 

(i)  From  coordinate  system  to  ’ ,  the  cone  is  tangent  to  the 
plane  yb ' ,  ofa ' ,  '  where  the  tangent  line  LpG  (see  Section3.2) 
is  coincident  with  y  '  axis.  Here,  we  have  (Fig.  3.8a). 


15 


cosa 


s  ma 


r- 


[M,  '  ]  =  -sina 
b  a 

0 

.  0 


cosa 

0 

0 


0  -dsina 

0  dtgasina 

1  0 

0  1 


(3.3.4) 


where  the  d  is  the  height  of  cone. 

(ii)  From  coordinate  system  S  '  to  S^'  the  tangent  line  of  cone 
and  plane  y ^  is  declined  with  an  angle  y  (see  Fig.  3.4),  also 
the  origin  o^'  and  o^  are  not  coincident.  Here  we  have  (Fig. 


3.8b). 


[M 


bb 


1 

0 

0 


L  o 


0  0 

cos  p  -siny 

siny  cosy 

0  0 


COSlp 


n 


COS  ip 


tgy 


(3.3.5) 


where  b  is  the  addendum  height  of  the  tool. 

(iii)  From  coordinate  system  to  Sa,  the  tool  is  declined  with 
a  pressure  angle  For  the  helical  gear,  the  i|»n  is  measured  at 
the  normal  cross  section.  Here,  we  have  (Fig.  3.8c). 


’cOSlJ)n 

-sim() 

n 

0 

0  ' 

Wab' 

= 

siml) 

n 

COSll> 

n 

0 

0 

(3.3.6) 

0 

0 

1 

0 

.  0 

0 

0 

1  J 

(iv)  From 

coordinate 

system 

^a 

to 

V 

the 

helix  angle  8p 

considered. 

Here  we  have  (Fig 

.  3. 

8d  ) 

rcos  8 

p 

0 

-sin6p 

0  ■ 

[H  ] 
ca 

= 

0 

1 

0 

0 

(3.3.7) 

sin8 

P 

cos  8 

P 

0 

i 

L  0  0  0  1  J 


[M  ] 


(3.3.7) 


I 


In  coordinate  system  Sc/  the  tool  surface  and  its  normal  can 


be  represented  as: 


"  ^a^b^b’^bd^d1 


(3.3.8) 


lnc  1  =  [LcallLab^Lbb-1[Lbd]tnd] 


where  3  x  3  L  matrix  is  from  corresponding  4  x  4  M  matrix, 


excluding  the  last  row  and  last  column.  Substituting  Eg. 
(3.3.1),  (3.3.3),  (3.3.4),  (3.3.5),  (3.3.6)  and  (3.3.7)  into 
Eq.  (3.3.8),  finally  we  obtain  the  tool  surface  and  its  normal 
in  coordinate  system  Sr  as 


where 


[n(p)] 

c 


!>>' 

!  cx 


cy 

n(p) 
L  cz  J 


(3.3.9) 


X(P)  =  u  cose  sina(cosacosi|j  cosB  +  sinasinij)  cosocosB 
c  P  P  n  p  n  p 

+  sinasinysinB  )  +u  sine  sina(sinj>  sinycosB  -  cosysinB  ) 

P  P  P  n  p  p 

+  u  cosa  (cosas in*  cosycosB  -sinacosti  cosB  +  cosasinysinB  ) 
p  n  p  rn  p  p 

-  (— —  -  - - — ■ — )  [sin<J»  cos y cos B  +  sinysinB  ] 

cosa  cosiji^cosy  n  p  p 


=  u  cose  sina  (cosasinib  -  sinacosib  cosy) 
p  p  n  rn 

-  u  sine  sinacosip  siny  -  ucosa  (sinasimb 
p  p  n  n 

,  d  b 

+  ( -  -  - ; - )  COSy  COS  lb 

cosa  cost^cosy  yn 


cos  acos  ipncos  y  ) 
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Z^^=u  cos0  sinct  (cosacosip  sing  +sinasintl»  cosysinR  -sinasinycos 6  ) 
cpp  n  p  *n  p  p 

+u  sin0  sina  (simji  sinysing  +  cosycosg  ) 
p  p  n  p  p 

+u  cosa(cosasim|i  cosysing  -  sinacosij)  sing  -  cosa  siny) 


+  ( 


cosa 


, — — — — )  (-cosysing  sin^  +  sinycosg  ) 
cosi|i  cosy  p  n  p 


(P)  _ 


cx 


cos0  cosa  (cosacostb  cosg  +  sinasintjt  cosycosB 
p  n  p  n  o 


+  sinasiny sing  )  +  sine  cosa(sintl)  sinycosg  -  cosysing  ) 
P  P  n  M  p  p 

+  sina  ( sinacos^cos g^-cosas  im^cosy  cos g^-cos asinysinfi^) 


=  cos  0  cosa  (cosas  intli  -sinacosij)  cosy)-sine  cosacosij)  sing 
cy  p  Tn  vn  p  n  p 

+  sina  (sinasinit  +  cosacosd)  cosy) 
n  Tn 


(P>  _ 


cz 


cos 0  cosa  (cosacostj/  sing  +sinas inti  cos ys  in g  -sin as  in y cos  g  ) 

ir  *  r  **  r  r 

+  sine  cosa(sint|j  sinysing  +  cosycosg  ) 
p  n  p  p 

+  sina  (sinacosij)  sing  -cosasinij)  cosysing  +cos  as  in  ycos  g  ) 
n  p  n  p  p 


3.4  Pinion  Tooth  Surface 

The  equation  of  meshing  of  the  generating  surface  and  the 
helical  pinion  tooth  surface  is  represented  by 


N(P)  .  V(PD  .  N(P).  (V(P)  _  v(l)}  .  o 


(3.4.1) 


We  can  also  use  equation  based  on  the  fact  that  the  contact 
normal  of  generating  and  generated  surfaces  must  intersect  the 
instantaneous  axis  of  rotation  I-I  (Litvin,  1968).  Thus,  we 
obtain  (Fig.  2.2). 
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(3.4.2) 


where  (Xc,  Yc,  ZQ)  are  the  coordinates  of  a  point  that  lies  on 
axis  I-I;  xc  ,  y  ,  are  the  coordinates  of  cone  surface; 

ncx'  ncy  an<^  ncz  are  tlie  projections  of  the  surface  unit 
normal.  From  Fig.  3.5,  it  is  known  that 


Equation  (3.4.2),  (3.4.3),  and  (3.3.9)  yield 


f,(u  ,9  , )=-u  [cose  (cosycosB  +sinysin\ii  sinB  ) 

Ipppp  p  p  np 

+sine  (sinasinycosB  -cosacosa>  sing  -sinasintfj  cosysinB  )] 
P  P  n  p  n  p 

+  "  ^oTilTcosy0  I^osepcos2a+sin2a).  [cosycosBp 

+sini|>  sinysinfi  1 -sine  cosacosi[<  sinB  ] 
n  p  p  n  p 

+r<)>  [cose  cosa (cosasimj/  -sinacos^  cosy) 
p  p  n  n 

+sina  (sinasinip  +cosacositi  cosy)-sin6  cosacosih  sintb  ] 
n  n  p  n  n 

-  0 


(3.4.4) 


The  helical  pinion  tooth  surface  can  be  obtained  by  transforming 
the  generating  tool  surface  E  to  coordinate  system  Sj,  together 
with  equation  of  meshing.  The  coordinate  transformation  in 
transition  from  Sc  to  Sf  is  represented  by  the  matrix  Mfc  as 
(Fig.  3.5) 

1  0  0  -r .  $ 

1  P 

[M  1  =0  1  0  0  (3.4.5) 

0  0  10 
.0  0  0  1 
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The  coordinate  transformation  in  transition  from  Sf  to  is 


represented  by  the  matrix 

as 

(Fig. 

3.5). 

COS#p 

sin#„ 

P 

0 

0  ~ 

“hf1  = 

-sin# 

P 

COS# 

P 

0 

0 

(3.4.6) 

0 

0 

1 

0 

1 

0 

0 

0 

1  . 

Therefore,  the  helical  pinion  surface  can  be  represented  in 
coordinate  system  as: 


lrl'  -  lMl=IIrcP)|  ■  i”1£ltM£cHr'P)] 


(3.4.7) 


together  with  Eq.  (3.4.4).  Substituting  Eq.  (3.4.5),  (3.4.6)  and 
(3.3.9)  into  Eq.  (3.4.7)  we  have 


[til 


( P ) 

*p)cos*p+yc  sin<(,p 
( P) 

iVsinVyc  cos4,p 


i 


(  3.4.8) 


Eq.  (3.4.8)  and  (3.4.4)  represents  the  pinion  tooth  surface 

where  x^  ,  y^  and  zc  are  expressed  m  Eq.  (3.3.9). 

Using  the  same  approach,  the  unit  normal  pinion  tooth 

surface  can  be  represented  by  Eq.  (3.4.9)  and  (3.4.4) 

(  P )  ( P )  (  p ) 

where  n  ,  n  and  n  are  expressed  in  Eq.  (3.3.9) 

c  x  cy  c  z 


[n1] 


n'  cos#„  +  n'  sin<fc 
cx  Tp  cy  p 

(P)  •  (P) 

ncx  Sln*p  +  ncy  cos*p 

«(P) 


(3.4.9) 
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It  is  clear  that  if  we  set  =  0,  the  pinion  becomes  spur 
pinion.  Therefore,  crowning  spur  pinion  using  conic  tool  is  only 
a  special  case  of  crowning  helical  pinion. 


3.5  Condition  of  Pinion  Non-Undercutting 

The  problem  of  undercutting  of  the  helical  pinion  tooth 
surface  by  crowning  is  related  with  the  appearance  on  the  pinion 
tooth  surface  of  singular  points.  From  differential  geometry,  it 
is  known  that  the  surface  point  is  singular  if  the  surface  normal 
is  equal  to  zero  at  such  a  point. 

Litvin,  proposed  a  method  to  determine  a  line  on  the  tool 
surface  whose  points  will  generate  singular  points  on  the  surface 
generated  by  the  tool.  This  line  designated  by  L  (Fig.  3.9)  must 
be  out  of  the  working  part  of  the  tool  surface  to  avoid 
undercutting  of  the  pinion  by  crowning. 

The  limiting  line  L  of  the  tool  surface  is  determined  by  the 
following  equations 


(w 


(3.5.1) 


fl(Up' 


(3.5.2) 


F(up, 


(3.5.3) 


Vector  equation  (3.5.1)  represents  the  tool  surface  [see  Eq. 
(3.3.9)];  equation  (3.5.2)  is  the  equation  of  meshing  [see  Eq. 
(3.4.4)]  and  equation  (3.5.3)  comes  from  the  requirement  of 
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limiting  line  L  as: 
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where  is  the  relative  velocity  of  the  tool  and  generated 

surface  at  generating  point,  represented  if  coordinate  system 

(cl) 

S~.  Actually,  we  can  write  V  as 


v<cl> 


(3.5.5) 


where  V^c ^  is  the  velocity  of  the  cutter  and 
of  the  pinion.  From  Fig.  2.2,  we  get 


ds_~ 

dt 

“p 

= 

0 

= 

0 

0 

.  0 

is  the  velocity 


(3.  5.6) 


v*1*  = 

U) 

x  r  ^  +  0  0 , 

X  CD 

~c 

~P 

~c  cl 

~P 

where 

id  = 

o  o 

o)  (Ts;  = 

rl^P' 
-r . 

~P 

.  1  . 

c  1 

1 

0 
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2  2  2  .  . 

X  =  [  ( AFsin  ct  -  CE)cos9p  +  ( ADs  m  ot  -  AEcos  a)sin9p 

+  ( AFcos 2a  -  DC)]  { Bcos0  +  Asin9  +  I) 

cr  IT 

Y  =  Dtgacos9  -  Ftgasin9  -  Ectga 

P  P 

Z  =  COSpCOSl^ 


23 


Up^9p^  >  cosa 


(3.5.11) 


The  analysis  of  Eq.  (3.5.10)  indicates  that  the  inequality 
(3.5.11)  is  satisfied  if  the  condition  of  helical  pinion  non- 
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undercutting  by  a  regular  rack  cutter  is  satisfied. 


3.6  Principal  Directions  and  Curvatures  of  Tooth  Surface 

A  simplified  approach  to  determine  principal  directions  and 

curvatures  of  helical  pinion  has  been  proposed  by  Litvin  (Litvin, 

1968).  The  main  idea  is  representing  the  principal  directions 

and  curvatures  of  generated  surface  E^  by  the  principal 

directions  and  curvatures  of  generating  surface  E^. 

Let  us  determine  the  principal  directions  of  curvatures  of 

tool  surface  e  .  The  tool  surface  is  a  cone  surface  and  its 
P 

principal  directions  coincide  with  the  direction  of  the  cone 
generatrix  and  the  direction  that  is  perpendicular  to  the  cone 
generatrix . 

The  Rodrigue's  formula  (Eq.  3.6.1)  can  be  used  to  obtain  the 
principal  curvature  and  directions 


< 


I, II 


V  = 
~r 


-n 

~r 


(3.6.1) 


where  V  is  the  velocity  of  a  point  that  moves  over  a  surface 

and  n  is  the  derivative  of  the  surface  unit  normal  n, 

~r 

when  n  changes  its  direction  due  to  the  motion  over  the  surface. 

Using  Eq.  3.6.1,  the  principal  directions  and  curvatures  of 
cone  surface  eq  can  be  expressed  in  coordinate  system  as: 
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(P) 


Si 


,(P)  - 

ii 


3~d  . 

~a 

'-sine  * 
0 

cose. 

(P) 

l 

30 

P 

30 

1  P 

r  K 

i 

u  cos  a 
P 

9rH 

~d  . 

3r  . 
~d 

'  cosesina ' 

• 

■ 

(p)  -  n 
<11  -  0 

3u 

P 

3U 

P 

“COSa 

.sinesina  . 

1' 

(3.6.2) 


The  negative  sign  for  indicates  that  the  curvature  center  is 
located  on  the  negative  direction  of  the  surface  normal.  The 
principal  curvatures  are  invariants  with  respect  to  the  used 
coordinate  system.  Whereas  the  principal  direction  will  be 
represented  in  coordinate  system  Sf  as 


(P 
I, II 


f 


lhc1ILcd>‘^}i]d 


(3.6.3) 


where  [ Lcd  1  =  [ Lca ^Lab^Lab '  ^Lb  ’d 1  '  3  x  3  L  matrices  can  be 
obtained  from  corresponding  4  x  4  M  matrices  [see  Eq.  (3.3.4), 
(3.3.5),  (3.3.6)  and  (3.3.7)],  L^c  is  3  x  3  unitary  matrix  (see 
Fig.  3.5). 

The  determination  of  principal  curvatures  and  directions  for 
the  pinion  tooth  is  based  on  the  following  equations  (see  Litvin, 
1968) 


tan  2a 


(PD  _ 


2b13b23 


h2  h2  /  (P> 

b23**b13~  KI 


-K(P)  )b 
KII  ,D33 


(3.6.4) 


( 1 ) _ <  1 )  _ 


II 


•ic 


H2  K  2 

b  23~b 13 


-  ((C(P)_K(P) 
y  I  *II 
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b  ^  3cos  2a 


(pl) 


(3.6.5) 
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(3.6.6) 
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(1) 

II 


+  K 


(P) 

II 


2 

23 


b 


33 


where  and  e!P\  el?^  are  the  principal  curvatures  and 

unit  vector  of  principal  directions  of  E^.  kii^  an(l 

e^1^,  are  the  principal  curvatures  and  unit  vectors  of 

principal  directions  of  Eq.  Angle  is  measured  counter 

clockwise  from  e  ^  ^  to  e^^  (Pig.  3.10).  The  coefficient  b^-j, 
b23  and  b33  have  been  derived  as  shown  in  (Litvin  1968)  but 
modified  for  the  case  when  a  rack  cutter  generates  a  gear.  The 
expressions  for  b-^,  b23»  and  b33  are  as  follows: 


b  13 

= 

’  (p) 

eII  *  ~p 

r  -K<p> 

ki 

0 

r  v(pi} 

1 

_ b  23 

-e'P>  .  « 

L  ~i  ~P  J 

0 

-K(p) 
KII  J 

V(PD 

<v 

e(p) 

*  ~I  J 

(3.6.7) 
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r  . .  ( p ) .  .  ,  „(pl),  (p )  .  (Dl )  (p).2 

=  [n  a)  V  c  ]  +  (n  to  V  ]  -  kt  (V  .  eT  ) 


/V*  (V 


p  ~tr 


-  4P)(V(P1)e(P))2  .  _ 
11  ~  ~  1  1 

m 


lp 


"ip  in  £<plV 


where  m 


lp 


d  (}> 
ds 


dm 


m 


IP 


l£ 


ds 


(3.6.8) 


The  vectors  used  in  Eq.  (3.6.7)  and  (3.6.8)  are  represent  in 
coordinate  system  Sf  (Fig.  3.5)  with  i ^ ,  i^ ,  k^  as  its  unit 
vectors  of  axes.  The  expressions  of  the  vectors  in  Eq.  (3.6.7) 
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and  (3.6.8)  are  as  follows  (Fig.  3.5) 


w  =  tu  kc  (3.6.9) 

~P  ~f 

is  the  angular  velocity  of  the  pinion  being  in  meshing  with  the 
rack  cutter. 


V(P) 

tr 


'wrlif 


(3. 6. 10) 


is  the  transfer  velocity  of  a  point  on  the  rack  cutter  that 
performs  translational  motion,  r^  is  the  radius  of  pinion 
centrode . 


r-v 

1 

o 

+ 

y- > 

1 

V<1>  =  -  X  r<?>  -  .  |  xf 

~tr  ~p  ~f  t 

=  0) 

x  -r,  <p 
c  1  p 

L  o  . 

0 

(3.6. 11) 


is  the  transfer  velocity  of  the  pinion 


v(pD  =  V(P)  _ 


-x  +r, 
c  Id 


0 


(3.6.12) 


is  the  "sliding"  velocity  -  the  velocity  of  a  point  of  rack 
cutter  with  respect  to  the  same  point  of  pinion. 

Substituting  Eg.  (3.6.9),  (3.6.10),  (3.6.11),  (3.6.12), 
(3.6.2)  and  (3.3.9)  into  Eq.  (3.6.7)  and  (3.6.8),  then 
substituting  Eq.  (3.6.7)  and  (3.6.8)  into  (3.6.4),  (3.6.5)  and 


(3.6.6),  finally  we  can  obtain  the  principal  directions  and 
curvatures  from  Eq.  (3.6.4),  (3.6.5),  (3.6.6)  and  (3.6.2).  Since 
the  expressions  of  tool  surface  and  its  principal  direction  are 
complicated  and  tedious.  It  is  difficult  for  us  to  write  down 
the  expression  of  e^  and  ej2^.  However,  a 

corresponding  computer  program  is  developed  for  determining  the 
principal  directions  and  curvatures  of  any  point  of  the  pinion 
surface  using  the  algorithm  discussed  above. 

The  same  approach  can  be  used  to  determine  the  principal 
directions  and  curvatures  of  regular  helical  gear.  However,  in 
this  case,  the  generating  surface  is  a  plane,  and  the  problem 
becomes  simpler. 

3.7  Contact  Ellipse  and  Bearing  Contact 

The  dimensions  and  orientation  of  the  instantaneous  contact 

ellipse  can  be  determined  based  on  the  equations  proposed  by 

Litvin.  (Litvin  1968).  The  input  data  for  computation  is: 

k ^ 1  ^ T  (i  =  1,2),  e!1*,  0(12)  and  e,  where  are  the  principal 

curvature  of  pinion  tooth  surface  E ^  and  gear  tooth 

surface  ei^  is  the  unit  vector  of  the  first  principal 

(  12 ) 

direction  of  jj.  a  is  the  angle  formed  by  the  unit  vectors 

of  principal  directions  e^^  and  e^2^  (Fig.  3.11)  and  e  is  the 
elastic  approach  of  the  contacting  surfaces.  The  bearing  contact 
is  formed  by  a  sot  of  instantaneous  contact  ellipses  that  move 
over  the  gear  tooth  surface  in  the  process  of  meshing. 

The  axes  of  the  contact  ellipse  are  directed  along 

the  n  axis  and  5  axis,  respectively.  The  orientation  of  contact 
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ellipse  is  determined  by  the  angle  a,  which  is  angle 
between  n  and  e^^.  (Fig.  3.11)  The  dimensions  of  ellipse  are 
determined  by  a*  and  b*  which  are  the  half  lengths  of  major  and 
minor  axis  of  the  ellipse  respectively.  The  following  equations 
are  used  to  determine  a,  a*  and  b*. 


A  = 

B  * 


1 

4 

1 

4 


w  tan(2n)  = 


r  <1}  In 

l<t  ~  K  t  '  lgl 


[<(1)~  *  (2) 
e  e 


g  ^sin2a 


+  1 9  ! 
(12) 


g1-g2cos2a 


(12) 


where 


(3.7.1) 


<(1)  =  A" 

€  I 


+  K 


9i  = 
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-  x 


(2) 
II 

(1) 

II  ' 


+  ic 


(2) 

II 


K 


(2) 

II 


and 


a* 


e 

!/2  b*  = 

e 

A 

B 

1/2 


(3.7.2) 


3.8  Simulation  of  Meshing  and  Determination  of  Transmission 

Errors 

The  simulation  of  meshing  is  a  part  of  the  computer  aided 
tooth  contact  analysis  (TCA)  program.  The  simulation  of  meshing 
is  based  on  equations  that  provide  the  continuous  tangency  of 
contacting  surfaces. 

To  simulate  the  meshing  of  a  crowned  helical  pinion  and 
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regular  involute  helical  gear  with  misaligned  axis,  we  will  use 
the  following  coordinate  system,  as  shown  in  Fig.  3.11  and  Fig. 


3.12:  (i)  Sf  is  rigidly  connected  to  the  frame  (ii)  an 
auxiliary  coordinate  system  that  is  also  rigidly  connected  to 
the  frame  (iii)  S-^  and  S2  are  rigidly  connected  to  the  pinion  and 
the  gear  respectively.  The  relations  between  S  ^  and  Sf  and 
between  S2  and  are  shown  in  Fig.  3.12  and  expressed  by  ^  and 
Mh2  as 


M 


f  1 


COS<t>  ^ 
-sin<)>  ^ 
0 
0 


sin$  ^ 
cos  4)  ^ 
0 
0 


0 

0 

1 

0 


0  ' 
0 
0 

1 . 


M 


h  2 


cos<J>2  -sin<ji2  0  0 

sin<)>„  cos4>9  0  0 

0  0  10 
0  0  0  1. 


(3.8.1) 


(3.8.2) 


It  is  obvious  that  the  rotation  axis  of  the  pinion  and  gear  are 
Zf (Zx)  and  Z^(Z2)  respectively.  Therefore,  the  error  of  assembly 
of  gears  can  be  simulated  with  orientation  and  location  of 
coordinate  system  with  respect  to  S^.  Figure  3.13  shows  the 
orientation  of  and  S^.  When  the  pinion  and  gear  axes  are 
crossed  (Fig.  3.13a)  and  intersected  (Fig  3.13b)  with  operating 
center  distance  C  =  0^0^.  It  must  be  emphasized  that  C  can  be 
different  from  the  nominal  center  distance  given  by  C°  =  r^  +  r2, 
where  r^  and  r2  are  the  radii  of  the  pinion  and  gear  pitch 
circles.  The  coordinate  transformation  from  Sh  to  Sf  is 
represented  by  [Mfhl  as  follows  (Fig.  3.13) 
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-cos Ay 
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sinAy 

0 

cosAy 
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0  ' 
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(3.8.3) 


where  the  gear  axes 


crossed  as  shown  in  Fig.  3.13a. 


-1 

0 

0 

0 


0  0  0 

-cosa^  -sinAj  C 
-sinAf  cosAf  0 

0  0  1 


(3.8.4) 


where  the  gear  axes  are  intersected  as  shown  in  Fig.  3.13b. 

The  helical  pinion  tooth  surface  and  its  normal  can  be 
expressed  in  coordinate  system  as: 


]  -  [Mfl 

]  [r1] 

(3.8.5) 

4-1 

II 

Hn,] 

where  ^ 

is  3  x 

3  matrix  from  M^. 

The  helical  gear  tooth 

surface  and 

its  normal  can  be  expressed 

in  coordinate  system 

cis 

^2> 

]  =  [Mfh 

1  ^ r 2  ^ 

(3.8.6) 

(n«2» 

1  =  [Lfh 

] [Lh2^  tr25 

where  and  is  3  x  3  matrices  from  and  M^2» 

For  a  helical  gear  (lefthand),  [  r  2  3  ar>d  can  be 

represented  as: 
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The  contact  of  gear  tooth  surface  is  simulated  in  the  TCA  program 
by  the  following  equations 
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E£1)(VV*1> 


(3.8.9) 


2f  )(<^p'9p^l) 


~f  ^  ^G'fcG'*2* 


(3.8.10) 


Vector  equation  (3.8.9)  comes  from  the  equality  of  the  position 
vectors  at  the  contact  point  and  provide  three  independent 
equations  whereas  vector  equation  (3.8.10)  comes  from  the 
equality  of  surface  unit  normals  and  provide  only  two  independent 
equations.  Therefore,  Eq. (3,8.9)  and  (3.3.10)  yield  five 
independent  scalar  equations  as  follows: 


f  i  ( 41 1  /  9p,<^p,^2,^G,<^G^  -  ^  i  ~  1,2,. ..,5)  (3.8.11) 


Equation  system  3.8.11  is  the  expression  in  implicit  form  of 
functions  of  one  variable,  i.e.  ^ .  Using  theorem  of  implicit 
function,  we  can  obtain 


2  =  *2U1} 


(3.8.12) 


considering  at  the  neighborhood  of  P°  =  {()>..,  0^4>  ,  4>_,tp,  <t>„)  where 


1  p*-Tp  2  G  G 


P°  satisfy  Eq.(3.8.11)  and 


D^l'^2f^3'^4'^5^ 

D(<f,p,0p,^2,tG'^G)^ 


(3.8.13) 


The  function  of  transmission  error  can  be  determined  by 
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1 


(3.8.14) 


A 2 ^  ^  i ^  ~  ~  N~  ^ 


3.9  Modification  of  Generating  Surface 

Using  cone  as  a  tool  surface  to  generate  the  helical  pinion 

is  good  way  to  localize  bearing  contact.  By  using  the  TCA 

program,  it  can  be  shown  that  the  transmission  errors  caused  by 

gear  misalignment  are  on  a  very  low  level.  However,  the  shape  of 

the  function  of  transmission  errors  is  unfavorable  as  shown  in 

Fig.  3.14(a).  This  shape  of  the  function  of  transmission  errors 

will  result  in  interruption  and  interference  with  change  of 

meshing  gear  teeth.  A  more  favorable  shape  of  the  function  of 

transmission  errors  is  shown  in  Fig.  3.14(b).  To  obtain  this 

shape  of  transmission  error  function,  we  need  to  modify 

generating  surface  T.  into  a  revolute  surface.  The  reason  for 

P 

using  revolute  surface  is  similar  to  the  case  for  crowning  a  spur 
pinion  (see  Litvin  and  Zhang,  1987). 

The  surface  of  revolution  is  generated  by  an  arc  circle  with 
radius  p.  The  arc  has  a  common  normal  with  the  cone  generatrix 
line  at  the  point  M  (Fig.  3.15(a)).  The  circular  arc  and  its 
normal  are  expressed  in  an  auxiliary  coordinate  system  S0  as 
follows 


X  =  p  [cos  a  +  B)-cosal  +  (— — —  -  sina 

e  cosa  cos<(i  cosy 

n 


=  (—"—  -  ) s i n a  -  2psin  sin(a  +  4) 

cosa  cos\|>ncosu  2  2 


Y  =  p[sin(a+8)  -  sina]  +  —  — — ; — -  cosa 

e  1  cosi^cosy 


(3.9.1) 
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sin(a+B  ) 
0 


(3.9.2) 


The  surface  of  revolution  is  generated  by  rotation  of  the 
circular  arc  about  the  y0  axis  and  can  be  represented  in 
coordinate  system  as  follows 


[rj  =  tLaJ[rJ  '  =  tLaJ[rJ 

d  de  e  d  de  e 


(3.9.3) 


where  (Fig.  9.2(b)) 
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(3.9.4) 


’cosepCOS(a  +  8) 

[nj  =  sin(a  +  8) 
a 

sine  cos(a  +  8) 

P 

The  installment  of  the  generating  surface  in  coordinate 
system  is  the  same  as  described  in  Section  3.3.  The  generated 
crowning  pinion  tooth  surface  can  be  obtained  using  the  same 
approach  as  that  described  in  section  3.4. 

The  meshing  of  gears  using  the  crowning  method  described  in 
this  section  has  been  simulated  by  numerical  methods.  The 
results  of  the  investigation  are  illustrated  with  the  following 
example . 

Given:  number  of  pinion  teeth  =  20,  number  of  gear  teeth 
N 2  =  40,  diametral  pitch  in  normal  section  Pn  =  10  in--*-,  pressure 
angle  in  normal  section  iJj^  =  20°,  helical  angle  8  =  15°.  The 
pinion  tooth  is  crowned  by  revolute  surface  with  generatrix 
arc  p  =  30  in.  The  revolute  surface  is  deviated  from  a  cone 
(comparing  in  figure  3.2(b)  and  (c)).  The  cone  has  half  apex 
angle  a  =  20°  and  bottom  radius  R  =  10.6  in. 

The  topology  of  the  pinion  tooth  surface  provides  a 
parabolic  type  of  predesigned  transmission  errors  with  d  -  6  arc 
seconds  (Fig.  2.3  (a))  and  a  path  contact  that  is  directed  across 
the  tooth  surface  (Fig.  3.1). 

The  influence  of  gear  misalignment  has  been  investigated 
with  the  developed  computer  program  and  the  results  of 
computation  are  represented  in  table  3.9.1  and  3.9.2  for  crossed 
and  intersected  gear  axes,  respectively.  The  misalignment  of 
gear  axes  is  5  arc  minutes. 
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The  results  of  computation  show  that  the  resulting  function 
of  transmission  errors  is  a  parabolic  one.  Thus  the  linear 
function  of  transmission  errors  that  was  caused  by  gear 
misalignment  has  been  absorbed  by  the  predesigned  parabolic 
function.  Fig.  3.14  show  the  results  of  transmissions  errors  for 
crossed  helical  gears  with  (Table  3.9.1)  and  without  (Table 
3.4.1)  crowning  of  the  pinion. 


TABLE  3.9.1  -  TRANSMISSION  ERRORS  OF  CROSSED  HELICAL  GEARS 


♦l 

(deg) 

-14 

-11 

-8 

-5 

-2 

1 

4 

_ 

A*, 

(sec ) 

-4.99 

-1.51 

0.65 

1.51 

1.05 

-0.75 

-3.87 

TABLE  3.9 

.2  -  TRANSMISSION  ERRORS  OF 

INTERSECTED  HELICAL  GEARS 

Uec,) 

mm 

-8 

-5 

-2 

1 

4 

7 

(sec  ) 

-6.15 

-2.  72 

o 

• 

O 

1 

0.  20 

. 1 

1 

-0.32 

-21.9 

-5.40 

4.  CROWNED  HELICAL  PINION  WITH  LONGITUDINAL  PATH  CONTACT 
4.1  Basic  Concepts  and  Considerations 

A  longitudinal  path  of  contact  means  that  the  gear  tooth 
surfaces  are  in  contact  at  a  point  at  every  instant  and  the 
instantaneous  contact  ellipse  moves  along  but  not  across  the 
surface  (Fig.  4.1(a)).  It  can  be  expected  that  this  type  of 
contact  provides  improved  conditions  of  lubrication.  Until  now 
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only  the  Novikov-Wildhaber 1 s  gears  could  provide  a  longitudinal 
path  of  contact.  A  disadvantage  of  this  type  of  gearing  is  their 
sensitivity  to  the  change  of  the  center  distance  and  the  axes 
misalignment.  The  sensitivity  to  non-ideal  orientation  of  the 
meshing  gears  cause  a  higher  level  of  gear  noise  in  comparison 
with  regular  involute  helical  gears.  Litvin  et  al.  (Litvin, 
1985)  proposed  a  compromising  type  of  non-conf ormal  helical  gears 
that  may  be  placed  between  regular  helical  gears  and  Novikov- 
Wildhaber  helical  gears.  The  gears  of  the  proposed  gear  train 
are  the  combination  of  a  regular  involute  helical  gear  and  a 
specially  crowned  helical  pinion.  The  investigation  of 
transmission  errors  for  helical  gears  with  a  longitudinal  path  of 
contact  shows  that  their  good  bearing  contact  is  accompanied  with 
an  undesirable  increased  level  of  linear  transmission  errors. 
The  authors  propose  to  compensate  this  disadvantage  by  a 
predesigned  parabolic  function  of  transmission  errors,  that  will 
absorb  the  linear  function  of  transmission  errors  (see  section 
2.3).  The  two  following  methods  for  derivation  of  the  pinion 
tooth  surface  will  the  modified  topology  will  now  be  considered. 


4.2  Method  1 

Consider  that  two  rigidly  connected  generating 

surfaces,  z_  and  Z  ,  are  used  for  the  generation  of  the  gear  and 
G  p  y 

the  pinion,  respectively  (Fig.  4.1(b)).  Surface  z~  is  a  plane 

Cj 

and  represents  the  surface  of  a  regular  rack-cutter; 

surface  z  is  a  cylindrical  surface  whose  cross-section  is  a 
P 

circular  arc.  We  may  imagine  that  while  surfaces  and  z 

G  p 
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translate,  as  the  pinion  and  the  gear  rotate  about  their  axes. 
To  provide  the  predesigned  parabolic  function  of  transmission 
errors  it  is  necessary  to  observe  the  following  transmission 
functions  by  generation 

N 

=  r2  -  const,  =  r2  (^  -  2a*p)  =  f<*p)  (4.2.1) 

Here:  u>G  and  u>p  are  the  angular  velocities  of  pinion  and  gear  by 

cutting;  V  is  the  velocity  of  the  rack-cutter  in  translational 
motion;  and  are  the  gear  and  pinion  tooth  numbers;  4>p  is 

the  angle  of  rotation  of  the  pinion  by  cutting.  The  generated 
gears  will  be  in  point  contact  at  every  instant  and  transform 
rotation  with  the  function 


<t>  2  ( <f|  i ) 


0  < 


(4.2.2) 


This  function  relates  the  angles  of  rotation  of  the  pinion  and 
the  gear,  ^  and  <t>2,  respectively,  for  one  cycle  of  meshing. 
The  predesigned  function  of  transmission  errors  is 


A*2  = 


(4.  2.3) 


It  is  evident  that  after  differentiation  of  function  (4.2.2)  we 
obtain  that  the  gear  ratio  satisfies  equation  (4.2.1), 

if  <ji^  and  <p  ^  are  used  instead  of  <j>p  and  $G. 

To  apply  this  method  of  generation  in  practice  it  is 
necessary  to  vary  the  angular  velocity  of  the  pinion  in  the 


40 


I 


process  of  its  generation  that  may  be  accomplished  by  a  computer 

controlled  machine  for  cutting. 

It  is  obvious  that  in  this  method  the  gear  tooth  surface  is 

kept  as  regular  skew  involute  surface  as  shown  in  Eg. (2. 4.1)  and 

(2.4.2)  since  its  generating  surface  l  is  a  plane  and  generating 

motion  is  V  =  <u  r  .  Now,  let  us  derive  the  pinion  tooth 

G  2 

surface  equations. 

The  pinion  generating  surface  E^  is  a  cylindrical  surface 
with  circular  arc  as  its  cross  section  and  can  be  represented  in 
an  auxiliary  coordinate  system  Sg  as  follows:  (see  Fig.  4.2) 


£a  = 


• 

X 

a 

R  [cos  (  ^n+a  )  -  cosij^] 

ya 

= 

R  [sin ( 4»n  +a  )  -  sintf/^] 

za 

t 

P 

.  1 

1 

(4.2.4) 


The  unit  surface  normal  is  represented  by 


cos ( ^n+«) 
sin(\|)n+« ) 


"  n 

r 

ax 

n  = 

n 

zz 

ay 

N 

(0 

c 

(4.2.5) 


From  coordinate  system  Sa  to  Sc  the  helical  angle  8p  is 
introduced.  Based  on  Fig.  3.8(d)  and  [McaJ  as  shown  in  Eq. 
(3.3.7),  we  have 


r(P>  = 
~c 


r 

c 

v(p) 

R[cos(t|in+a)  -  cosily]  -  tosin6p 

yc 

2  (P) 
c 

=  [Mca]£a  = 

R[sin ( \J)n+a )  -  sin i{/n ] 

R  [cos  (il«n+a)  -  cosi|»n]  +  tpSingp 

.  1 

1 

( 

(4.2.6) 
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J-  f 


(p>  _ 


[L  ]  r  = 
ca  ~a 


cos  ( ij)n+a  )cos  6p 

sin ( ^n+a ) 

cos  ( il>  +a)sinB 
n  p 


(4.2.7) 


The  generating  process  is  described  in  the  Fig.  3.5  where 


s  ■  r2»P<iq  -  aV 


(4.2.3) 


The  equation  of  meshing  of  the  generating  surface  and  the 
helical  pinion  tooth  surface  is  given  by: 


n(p)  .  v(pl)  =  n(p)  .  (V(P)  -  V(1))  =  0 


(4.2.9) 


In  the  system  shown  in  Fig.  3.5,  it  can  be  found  that 


V(P)  «  0 


-r.-y"” 

1  c 

r(l)  _  .  y(p) 

'  =  -S  +  X_  u). 


(4.2.10) 


Taking  into  account  that  -rr-2  =  u>  and  from  Eq.(4.2.1)  we  have 
*  dt  p 


+  r,+yip) 


-r 2 ( '  2aV  +  rl+V 


r(pD  =  V(P)  _  V(D  = 


-xcp)+r2*p<5r  -  "V 
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{ n  ) 

y/  +  2r_a< (> 

:c  2  Yp 


-x<P}  +  tr2cf>  (r— =• 
C  2yp  N 2 


a<b  ) 


(4.2.11) 


Substituting  Eq. (4.2.6),  (4.2.7)  and  (4.2.11)  into  Eq.  (4.2.9). 

the  equation  of  meshing  can  be  obtained  as: 


cos  (ii>n+a ) cos 8^ {  R[sin(ij>n+a)-sintlJn] 

+  sin(<J>n+a){  -R  [cos  (  )  -cosi^]  - 

=  0 


+  2r  a<f>  } 

2  p 

N1 

tsin$p  }  +  -  a*p) 

(4.2.12) 


Equation  of  meshing  gives  the  relation  among  three  variables, 

that  is,  t  ,  a  and  4>  It  is  obvious  that  from  Eq.  4.2.12,  the 
”  P 

equation  of  meshing  can  be  rewritten  as: 


tp  ~  cotan(^n+a)cotan(/?p)  {  R  [s  in  ( ^n+a  ) -sini)>n  ] +2r2a<j>p } 


+  Tikr  {-R[c°sUn+a)-cos*n]  +  r2*  -  a*  )} 

p  r  2  c 


(4.2. 13) 

The  pinion  tooth  surface  can  be  determined  with  the  same  approach 
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as  that  described  in  the  section  3.4  and  represented  by 
Eq. (3.4.8)  and  (3.4.9)  together  with  Eq. ( 4. 2. 6) , ( 4. 2. 7)  and 
(4.2.13). 

4.3  Method  2 

The  derivation  of  the  crowned  pinion  tooth  surface  is  based 
on  two  stages  of  synthesis.  On  the  first  stage  it  is  assumed  the 
pinion  tooth  surface  z ^  is  exact  conjugate  surface  to  the  gear 
tooth  surface  which  is  regular  skew  involute  surface  under  the 
condition  that  the  rotation  transformed  by  the  pinion  and  the 
year  is  described  in  Eq. (4.2.2)  with  predesigned  parabolic 
function  of  transmission  errors. 

On  the  second  stage  of  synthesis  it  is  necessary  to  localize 
the  bearing  contact  and  substitute  the  instantaneous  line  contact 
by  a  point  contact.  This  becomes  possible  if  the  pinion  tooth 
surface  will  be  deviated  as  it  is  shown  in  Fig.  4.1(c).  Only  a 
narrow  strip,  L,  will  be  kept  while  z  ^  will  be  changed  into  £ | . 
The  deviation  of  Z  j  with  respect  to  Z  ^  may  be  accomplished  in 
various  ways,  for  instance,  in  such  a  way,  that  the  cross-section 
of  is  only  a  circular  arc.  The  generation  of  rj  requires  a 
computer  controlled  machine  to  relate  the  motions  of  the  tool 
surface  and  the  being  generated  pinion  surface  zj.  The  tool 
surface  (it  may  be  only  a  plane)  and  z  |  will  be  in  point  contact 
in  the  process  of  generation. 

Now,  let  us  derive  the  pinion  tooth  surface  equations  based 
on  the  two  stages  discussed  above.  First,  we  can  derive  the 
intermediate  pinion  tooth  surface  £  ^  as  that  generated  by  gear 
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tooth  surface  r^.  The  relation  of  rotation  for  the  pinion  and 
the  gear  is  shown  in  Eq.(4.2.2).  The  gear  tooth  surface  an^ 
its  unit  normal  are  shown  in  Eq. (2.4.3)  and  (2.4.4). 


To  represent  the  gear  tooth  surface  1 2  in  the  fixed 
coordinate  system  Sf  in  Fig.  4.3,  it  is  necessary  to  transform 


(t- 2 ]  and  [n2] 

into  (rg2) 

1  and  [n<2)] 

as  : 

lr‘2,l  . 

=  lMf 21 lr2J 

tn'2)]  ■ 

=  tLf2l ^n2^ 

(4.3.1) 

I 

-cos  <j>  2 

s  i  n<f>  2  0 

O' 

where  [M^l  = 

-sin<(>  2 

-cos  <j>  2  0 

0 

0 

0  1 

C 

0 

0  0 

1. 

Substituting  Eq. (2.4.3)  and(2.4.4)  into  Eq. (4.3.1)  we  obtain: 

T  2 

COS 

cosUg-*q-42)  [-tGsinBp+r2»G)-r2sinUG-*2) 


[r<2>] 


x(2) 

xf 

Zf 

1 


2 

Cos  ijt 


C  '  CoTli  C°s(*G“*C  -*2)l_tGsin8p+r2*G1‘r2sin(V 


2  2  2 

tGcosSp(  1+sin  ^pn+g  8p)-r24Gsin  tj/^tg  8^ 


(4.3.2) 


tn‘2>l 


rn<2>  1 

nfx 

(2) 

nfy 

n(2) 

n£z 

L 

costf  cos6„cos  ($r-40)+sini|>  sin(  4»^— 4» -> ) 
n  p  \3  £.  n  kj  £. 

-cosi|;ncosSpSin(((iG-4i2)+sini()ncos  ( $G-4>2) 


costji  sins 
n  p 


(4.3.3) 
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The  generating  process  is  shown  in  Fig.  4.3 
where  4>  ^  and  <j>  is  given  by  Eq.  (4.2.2).  The  equation  of  meshing 
of  the  generating  surface  T.  ^  and  generated  surface  E  ^  is 
described  as: 


H<21  .  V<21>  =  N(2V2>  -  V(1))  =  0 


(4.3.4) 


The  relative  velocity  v  can  be  expressed  as 


r(  21)  =  (2)  _  (1)  = 


where 


jo 2  x  ^°2°1  +  Ef )  “  “i  x  £f 


(4.3.5) 


(1)  _ 


0  (it,  ;  to 

1  1  ~ 


(2)  _ 


w2  = 


0  to,  [ 

1  1 


1  [n2  2a  ^  1  ^ 


°2°1  =  -C  ; 


(2) 

Ef  =  yf 


and  the  relation  between  and  u  ^  is  obtained  by  taking 
derivative  of  Eq.  (4.2.2).  By  rearranging  and  simplifying  Eq. 


(4.3.5),  we  obtain 


V  =  X 


1  ^!'T'+  l~2a<t>^]  +  0  to  i  [^—  -  2a 

1  0  2 


(4.3.6) 


Substituting  Eg. (4.3.6)  (4.3.2)  and  (4.3.3)  into  Eq. (4.3.4)  and 


,6 
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m 


*  ■  » 


1 


t. 


▼ 


where 


2, 

COS  ip 

X1  =  ^sT  f  cos(^1+^c-X(^1))[-tGsinBp+r2^G] 

+  r2sin  ( ({>  -x  ( <p  1 )  )  -  Csin^ 


2 

COS  ij) 


Yi  = 


1  COSl 


”  sin((f>1+i|»c-X((f>1)  [-tGsin6  +r2*G] 


-  r2cos  ( <j>  ^-X  ( <j>  ^  )  )  +  Gcosij)^ 


.  2 


Z,  =  t~cos8  (  1+sin  \Ji  +g  8  )  -  r0d>„  sin  tot 
lop  np  Zvj  n 


g  p 


COS  1 p  cos  8 

n,  =  - - - ^  cos  ( 4> , -X  (  $  ,  )  +i|i  ) 

lx  COSlJ)  vyl  'u  yc 


n 


iy 


cosij)  cose 

_ n _ 2 

cosij/ 

c 


sinf^-A  ( 4>  1  )+^c) 


n.  =  cosi p  sin8 
lz  n  p 


Equations  (4.3.10),  (4.3.8)  and  (4.2.2)  represent  pinion 
tooth  surface  z^  in  coordinate  system  S-^.  But  z^  is  only  an 
intermediate  surface.  Now  we  come  to  the  second  stage,  that  is, 
to  deviate  i ^  into  zj.  The  surface  z|  must  satisfy  such 
condition  that  when  it  is  in  mesh  with  gear  tooth 
surface  z 2  without  misalignment,  the  contact  path  will  be  in 
longitudinal  direction.  The  surface  zj  will  be  formed  in  two 
steps:  (i)  the  contact  path  curve  in  z^  is  chosen  and  kept, 
(ii)  the  profile  of  surface  Z^  in  transverse  section  is  replaced 
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by  a  smaller  circular  arc  attached  on  the  L  in  such  a  way  that 
the  original  normal  of  along  L  is  kept. 

To  choose  curve  L  on  Ij  we  could  define  that  the  contact 
path  on  the  gear  tooth  surface  is  in  the  middle  of  the  tooth, 
that  is 


r 


(4.3. 11) 


where  r  =  /x^  +  y 2  is  the  distance  of  the  surface  point  to  gear 
axis  in  transverse  section.  Substituting  Eq.  (2.4.3)  into  Eq. 
(4.3.10)  and  simplyfing  it,  we  can  obtain  the  relation  of  surface 
parameter  along  the  contact  path  as 


r  ,,<{>_  -  t„sin8  =  0  (4.3.12) 

2  G  G  p 

Applying  Eq.(4.3.12)  to  Eq.(4.3.10),  we  can  find  the  curve  L  on 
pinion  tooth  surface  E ^  as: 

x^  =  r  2s  in  ( <t>  ( <f>  ^  ) )  -  Csin<j>-^ 

yl  =  _r  2COS  ^  1~X  ^  1  ^  +  Ccos<*’i 

N1  2 

z1  =  r2ctg8p(-j^-  <|>  1  -  a  4>  x  +  \  ( $  1 )  )  (4.3.13) 

where  the  4>  ^  is  the  curve  parameter  and  X  ( ^ )  is  described  in 
Eg .  (4.3.8). 
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Now  we  should  attach  the  circular  arc  to  the  L  described  by 
Eq.  (4.3.13)  in  the  transverse  section  to  form  new  surface  e^. 
Also  the  tangent  of  the  arc  at  the  point  on  the  line  L  must  be 
perpendicular  to  the  normal  of  E  ^  described  in  Eq.  (4.3.10).  As 
shown  in  Fig.  4.4,  the  equation  of  new  surface  e^  are: 

x^  =  x^  +  r  [cos  (  u  +a  )  -cosy  ] 

y^  =  +  R[sin(u+a)-siny]  (4.3.14) 


z 


1 


z 


1 


where  m  =  4>^  -  A($^)+ij>  .  Equation  (4.3.14)  describes  the  new 
pinion  surface  E ' .  In  Eq.  (4.3.14)  when  a  =  0,  the  designed 
contact  path  L  is  obtained. 


4.4  Discussion  and  Example 

In  section  4.2  and  4.3,  two  methods  have  been  presented  with 
derivation  of  the  equations  of  pinion  tooth  surface.  Comparing 
the  two  methods  for  the  generation  of  the  pinion  tooth  surface, 
it  may  be  concluded  that  both  provide  a  localized  bearing 
contact,  a  longitudinal  path  of  contact  and  predesigned  parabolic 
function  of  transmission  errors.  The  difference  between  these 
methods  is  that  the  tool  and  pinion  tooth  surfaces  are  in  line 
contact  by  applying  the  first  method  for  generation  and  in  point 
contact  by  the  second  one.  The  disadvantage  of  both  methods  for 
crowning  of  the  pinion  is  that  the  transmission  errors  caused  by 
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gear  misalignment  are  large  and  it  is  necessary  to  envision  a 
high  level  of  the  predesigned  parabolic  function  for  the 


absorption  of  transmission  errors.  This  is  illustrated  with  the 
following  example  (the  algorithms  for  simulation  have  been 
discussed  in  Section  3.8). 

Given  (the  data  is  from  ref.  3):  pinion  tooth  number  = 
12,  gear  tooth  number  N2  =  94;  diametral  pitch  in  normal  section 
Pn  =  2  in-^-;  pressure  angle  in  normal  section  =  30°;  helical 
angle  0  =  15° . 

The  pinion  tooth  surface  is  a  crowned  surface  whose  cross- 
section  is  an  arc  of  a  circle  of  the  radius  0.  3584in.  The 
predesigned  parabolic  function  is  of  the  level  d  =  25  arc  seconds 
(Fig.  2.3(a)). 

Consider  now  that  the  axes  of  the  gear  and  the  pinion  are 
crossed  and  the  crossing  angle  is  3  arc  minutes.  The  computer 
program  for  the  simulation  of  meshing  provides  the  data  of 
transmission  errors  that  is  given  in  Table  4.1.  The  data  of 
Table  4.1  shows  that  the  resulting  function  of  transmission 
errors  is  a  parabolic  function.  Thus,  the  linear  function  of 
transmission  errors  caused  by  misalignment  of  gear  axes  has  been 
absorbed  by  the  predesigned  parabolic  function. 

Table  4.2  represents  the  transmission  errors  for  the  same 
helical  gears  for  the  case  when  the  gear  axes  are  intersected  and 
form  an  angle  of  3  arc  minutes.  The  resulting  function  of 
transmission  errors  is  again  a  parabolic  function  with  the  level 
d  =  26.2  arc  seconds.  The  relatively  high  level  of  transmission 
errors  is  the  price  that  must  be  paid  for  the  longitudinal  path 
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of  contact.  However,  the  proposed  topology  of  the  pinion  tooth 
surface  provides  a  reduction  of  the  level  of  gear  noise  since  the 
linear  function  of  transmission  errors  is  substituted  by  a 


parabolic  function. 


TABLE  4.1  TRANSMISSION  ERRORS  FOR  CROSSED  HELICAL  GEARS 


♦l 

(deg ) 

-23 

a 

-13 

-8 

-3 

B 

7 

m 

-17.94 

-3.  06 

■ 

8.23 

4.84 

-19. 37 

TABLE  4.2  TRANSMISSION  ERRORS  OF  INTERSECTED  HELICAL  GEARS 


<t>  l 

(deg) 

-20 

-15 

-10 

-5 

0 

5 

10 

-23. 06 

-8.50 

L_ 

0.15 

2.  96 

0.  00 

-3.  66 

-22.95 

5.  DEFORMATION  OF  HELICAL  GEAR  SHAFT 
5.1  Basic  Concepts  and  Considerations 

Deformation  of  gear  shafts  always  exists  when  the  gears  are 
used  to  transmit  power.  This  is  because  under  the  load  the 
force  applied  on  gear  tooth  surface  is  transferred  to  the  gear 
shaft  and  the  shaft  is  not  rigid  body.  It  can  be  proven  that  the 
deformation  of  gear  shaft  results  in  the  same  effects  as 
misalignment  induced  by  assembly.  The  misalignment  from  assembly 
could  be  reduced  to  as  little  as  possible.  But  the  shaft 
deformation  is  inevitable.  Fortunately,  all  the  transmission 
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errors  and  shift  of  bearing  contact  due  to  deformation  of  gear 
shaft  can  be  compensated  by  crowning  helical  pinion  tooth  surface 
with  the  methods  discussed  in  Chapters  3  and  4. 


5.2  Force  Applied  on  Gear  Shaft 

When  the  pinion  and  gear  mesh  a  force,  between  their  tooth 
surfaces,  is  applied  on  the  contact  point.  The  direction  of  the 
force  is  along  the  normal  of  the  contact  point.  For  the  case  of 
regular  helical  gears,  in  the  fixed  coordinate  system  Sf  as 
described  in  section  3.8,  the  force  direction  is  a  constant  since 
the  normal  of  the  contact  point  is  a  constant  (Litvin,  1968)  and 
can  be  expressed  as: 


n  = 
~c 


'  COSlJl 

cos  8 

n 

sin* 

n 

COS(l 

n 

sine 

(5.2.1) 


As  shown  in  Fig.  5.1,  the  force  applied  on  the  pinion  tooth 
surface  is: 


F  =  F  +  F  =  -F 
~  ~t  ~z 


cos*  cos8„ 
n  p 

s  inty 


cos*  sin8 
n  p 


(5.2.2) 


cosiK  cos8„ 
n  p 

0 

where  F  =  -F 

sinij> 

n 

N 

II 

1 

^3 

0 

0 

cosij)  sine 

L  n  pJ 

i 
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F.  is  called  as  transverse  force  since  it  is  applied  on  the 
~t 

transverse  section  of  the  pinion  and  the  gear  and  F^  is  called  as 

axial  force.  Transverse  force  F  is  always  in  tangency  with  the 

base  circle  in  tranverse  section.  Therefore,  if  the  torque 

transferred  by  the  gears  is  constant,  the  magnitude  of  transverse 

force  is  also  constant.  So  is  the  magnitude  of  axial  force  since 

it  has  certain  ratio  with  transverse  force. 

Both  transverse  force  and  axial  force  can  be  transferred  to 

the  axis  of  gear  shaft  with  resultant 

torque  (designated  by  F*  and  F*)  .  For  the  transverse  force,  the 

~t  ~z 

resultant  torque  is  the  torque  transferred  by  the  gear.  For  the 
axial  force,  the  resultant  torque  is  balanced  by  the  support 
bearings.  Assuming  the  force  is  applied  on  the  middle  section  of 
the  gear,  after  the  force  is  decomposed  and  transferred  to  the 
axis  of  gear  shaft  as  shown  in  Fig.  5.1  (b),  both  the  transverse 
force  and  axial  force  will  act  on  the  origin  of  coordinate  system 
Sf.  Actually,  both  forces  will  cause  deformation  of  the  shaft. 
But  since  the  axial  force  only  cause  very  small  tension  or 
compression  of  the  shaft,  it  can  be  neglected. 

5.3  Modelling  of  Shaft  Deformation 

As  shown  in  Fig.  5.2,  the  transverse  force  F  is  applied  on 
the  point  A  which  is  the  center  of  the  shaft  corresponding  to  the 
pinion  or  gear  middle  cross  section.  Under  the  force,  the 
deformation  of  the  shaft  is  composed  of  two  parts,  that  is, 
deflection  of  shaft  at  the  point  A  designated  by  Vp  and  rotatior 
of  the  shaft  cross  section  with  A  as  a  center  designated  by 
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Ap.  The  value  Vp  and  Ap  depend  on  the  magnitude  of  transverse 
force,  geometry  and  material  of  shaft  and  the  way  how  the  shaft 
is  supported.  (Timoshenko  1973). 

Now,  let  us  model  the  deformation  of  the  helical  pinion 
shaft.  For  convenience  and  consistency  with  the  previous 
chapters,  We  establish  our  coordinate  system  as  follows: 

(1)  As  shown  in  Fig.  5.3a  3^  is  a  fixed  coordinate  system  and 
rigidly  connected  to  the  frame.  Sg  is  an  auxiliary  coordinate 
system  with  the  Y_  axis  coincident  with  the  direction  of 

Cl 

transverse  force.  The  angle  is  the  pressure  angle  in  the 
transverse  section  of  helical  gear.  The  matrix  [Mfal  can 
transfer  vector  from  Sa  to  Sf  and  is  expressed  as: 


[Mfa]  = 


sini|i( 

-COSl)^ 

0 

0 


COSlfl( 

sin*J>c 

0 

0 


0 

0 

1 

0 


0 

0 

0 

1 


(5.3.1) 


(2)  Coordinate  system  Sfl  and  Sa,  [Fig.  5.3(b)]  are  connected  to 
the  axis  of  pinion  shaft.  Without  deformation,  '  and  Sa '  are 
coincident  with  Sf  and  Sa  respectively.  The  matrix  [Ma  *  f  *  1  can 
transfer  vector  from  '  to  Sg i  and  is  expressed  as: 


‘Vf1  - 


sinifi^ 

COSlj)( 

0 

0 


-COSip 

sin<|> 

0 

0 


0 

1 

0 


0 

0 

1 


(5.3.2) 


(3)  Coordinate  systems  Sa  and  Sa '  are  not  coincident  when  shaft 
deformation  occurs.  As  shown  in  Fig.  5.2,  their  relation  can  be 
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expressed  as: 


M 


aa  ' 


* 

1 

0 

0 

0 


0  0  0 

cosX  -sinX  -v 

P  P  P 

sinX  cosX  0 

P  P 

0  0  1 


(5.3.3) 


(4)  As  described  in  previous  chapter,  pinion  tooth  surface  can 
be  expressed  in  coordinate  systems  S-^  as  a  position 
vector  with  its  normal  n^.  The  relation  of  Sj  and  '  is 
shown  in  Fig.  5.3(c)  and  expressed  as: 


M 


f '1 


COSij)  ^ 

-s  i  n  <j>  ^ 
0 
0 


sin<j>i 
COSij)  ^ 
0 
0 


0  0  ‘ 

0  0 
1  0 
0  1 


(5.3.4) 


After  the  coordinate  systems  are  established,  it  is  easy  to  find 
that  the  deformation  of  the  pinion  shaft  can  be  modelled  by 
matrix  niff'  as: 


[Mff,l  =  (MfJ  [M_J  [MaIfI] 


fa 


aa 


a '  f 


1+cos  tp  (cosX  -1)  ,siml>  cosij)  (cosX  -1)  ,  -cosij)  sinX 

C  p  O  C  p  q 

2 

sinij;  cosij/  (cosX  -1)  ,  1+sin  if;  (cosX  -1)  ,  -sinij/  sinX 

c  c  p  c  p 


sinX  cosij) 

P  c 


sinX  s imb 
P  c 


c 

cos  X 


,  -vpcos^c 
,  -vpsin^c 
,  0 
1 


(5.3.5) 


56 


Since  X 

P 

of  sinX 

P 

written 


is  very  small  angle,  it 
and  one  instead  of  cosX  . 


as : 


P 


is  reasonable  to 
Therefore,  Eq. 


use  X 

P 

(5.3.5) 


instead 
can  be 


[M 


ff 


J 


1 

0 

X  COSlj)  , 

p  G 
0 


0 

1 

XpSin^, 


0 


-x  cosij) 

P  c 

-X  sinij>  , 
P  c 

1 

0 


-V  cos\l 
P  C 

-V  sintl; 
p  vc 

0 

1 

(5.3.6) 


Also,  for  the  gear  shaft,  the  coordinate  systems  SG,  SG,, 

S^,  and  S^'  are  used  instead  of  Sf,  Sf '  ,  Sa  and  Sa', 

and  X  _ ,  and  V_  are  used  instead  of  X  and  V  .  Then,  using  the 
G  G  p  p 

same  ideas,  the  deformation  of  gear  shaft  can  be  modelled  by 
matrix  [maa • ]  as: 


GG' 


l+cos^ib  (cosX„-l)  ,  sini|>  cosij>  (cosX„~l),  cosiji  sinX~  ,  V  cosiJj 

Co  CCvj  CvjVj 

2 

sin<|>  cosip  (cosX„-l),  1+sin  i|>  (cosX„-l)  ,  cosi|j  sinX„  ,  V  cosi|) 

CCo  Cvj  Coo 


-sinX„cos\|> 

G  c 


,  -sinXGsim|>c 
,  0 


cos  X, 


0 

1 


Or  considering  that  X  is  very  small  angle: 

O 


(5.3.7) 


^MGG' ^ 


1 

0 

XGcos*c 


0 

1 

-X_sini(i 

O  ( 

0 


+  XgCOSl(ic 

+XGsin^c 

1 

0 


VGCOS*c 

VGS in  ^c 
0 

1 


(5.3.8) 


It  must  be  emphasized  that  the  gear  tooth  surfaces  are  expressed 


57 


in  coordinate  system  S2.  The  relation  of  S2  and  Sf '  is  shown  i 
Fig.  5.4  and  expressed  as: 

-cos^  sin<|)2  0  0 

[Mq  ,  2}  =  -sin$2  -cos4>2  0  0  (  5.3.9) 

0  0  10 

0  0  0  1 


Also  as  shown  in  Fig.  5.5,  the  coordinate  system  Sf  and  SG  an 
not  coincident.  There  is  a  center  distance  C  between  thei: 
origins.  Therefore  M^G  is  represented  as: 


(5.3.10 


Now  it  is  easy  to  simulate  the  performance  of  the  gears  with  thi 
deformation  of  their  shafts.  Assuming,  there  is  a  helical  pinio 
in  coordinate  system  represented  by  position  vector  r^  an< 
normal  vector  n^  and  a  helical  gear  in  coordinate  system  S 
represented  by  and  n2»  we  can  write  the  tooth  contact  equatio 
as : 


where  L  matrices  are  3x3  matrices  from  corresponding  4x4 
matrices,  crossing  4th  row  and  4th  column  and  M  matrices  can  t 
found  in  this  section.  Applying  the  same  ideas  discussed  i 
Chapter  2  and  used  in  section  3.8,  we  can  find  the  transmissic 
errors  and  the  shift  of  the  bearing  contact  by  computer  aide 
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simulation. 


It  is  interesting  to  mention  that  the  deformation  of  the 

shaft  can  be  expressed  as  a  combination  of  misalignment  with 

crossing  axes,  misalignment  with  intersecting  axes,  and  change  of 

the  center  distance.  This  is  because,  for  example 

1  0  -A  COS! I)  -V  cos * 

pc  p  c 

[MCCI]=  0  1  -A  sin*  -V  sin* 

ff  p  c  p  yc 

A  cos*  A  sin*  1  0 

p  c  p  c 

0  0  0  1 


1  0  -A  cosiJj  0  1  f  1 

P  C 

0  10  0  0 
A  cos\(i  0  1  0  0 

P  0  °  0  0 


0  0 

1  -A  sin* 

p  yc 

A  simj;  1 

P  c 

o  o 


10  0 
0  10 
0  0  1 
0  0  0 


-V  cos* 
p  c 

-V  siniji 
P  c 

0 

1 


(5.3.12) 


where  the  three  decomposed  matrices  are  represented  by  the  matrix 

for  crossing  axis  with  small  angle  A^cosi|>c,  matrix  for 

intersecting  axis  with  small  angle  A  sinij;  ,  and  matrix  for  axis 

pc 

displacement . 


5.4  Example  and  Discussion 

The  results  of  investigation  in  this  chapter  are  illustrated 
with  the  following  example.  Given:  number  of  pinion  tooth 
N}=  20,  number  of  gear  tooth  N2  =  40,  diametral  pitch  in  normal 
section  Pn  =  10  in~l,  pressure  angle  in  normal 
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section  ip  =  20°,  helical  angle  8  =  15°,  gear  tooth  length  L  = 
5/Pn.  Also  assume  deformation  values  Vp  »  Vq  =  0.0125  in, 
ip  =  =  2  min.  The  computer  program  for  simulation  provides 

the  data  of  transmission  errors  in  Table  5.1 

From  Table  5.1,  it  is  known  that  the  transmission  errors  due 
to  gear  shaft  deformation  is  an  approximately  linear  function  for 
regular  skew  involute  pinion  and  gear.  Similar  to  the  case  of 
gear  axis  misalignment,  the  linear  function  of  transmission 
errors  can  be  absorbed  by  predesigned  parabolic  function  of 
transmission  errors  that  is  obtained  by  crowning  helical  pinion 
tooth  surface. 

TABLE  5.1  TRANSMISSION  ERRORS  FOR  HELICAL  GEAR  WITH 
DEFORMED  SHAFT 


(Aeg) 

3 

6 

9 

12 

15 

8 

21 

( sec  ) 

1.38 

3.  16 

4.74 

6.32 

7.90 

9.48 

11.06 

6.  CONCLUSION 

Several  methods  of  crowning  helical  pinion  tooth  surface 
have  been  developed.  The  modified  pinion  tooth  surface  cai 
provide  predesigned  parabolic  function  of  transmission  error: 
that  are  able  to  absorb  linear  function  of  transmission  error 
induced  by  misalignment.  Also,  the  modified  pinion  tooth  surfac 
can  improve  the  bearing  contact.  Principles  of  computer  aide 
simulation  of  meshing,  contact,  and  respective  computer  program 
have  also  been  developed.  The  numerical  results  of  examples  o 


cronwed  helical  pinion  in  mesh  with  regular  helical  gear  show 
that  the  ideas  of  crowning  are  useful  to  get  favourable  bearing 
contact  and  allowable  transmission  errors.  But  the  synthesis  of 
pinion  tooth  surface  should  be  based  on  a  compromise  between  the 
requirements  of  transmission  errors  and  the  patterns  of  the 
bearing  contact. 
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FIG.  3.2  Generating  Surfaces 


FIG. 


FIG.  3.4  Generating  Tool  Surfaces 
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3.6  Pinion  and  Gear  Generatina  Surfaces 
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FIG.  3.12  Coordinate  Syst 
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FIG.  3. 
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13  Coordinate  Systems  for  Simulation  of  Gear  Misalignme 
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FIG.  4.2  Helical  Pinion  Generating  Surface 
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Ya 


FIG.  5.2  Shaft  Deflection 
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FLOWCHART  FOR  PROGRAM  I 


% 


Q  START  ) 


/  INPUT  GIVEN  DATA  7 


1  FIND  AUXILIARY  VALUES 
|  FOR  CALCULATION 

CALCULATE  PROFILE  OF 
PINION  SURFACE  CUT 

BY  PLANE  Z  =  CONST. 

COMPARE  THE  PROFILE 
WITH  INVOLUTE  CURVE 
AT  MIDDLE  SECTION 

_ 

(  STOP  ) 
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L 


<p.  . . # . 


o  o  o  n  n 


******************************************************* 


i 


PROGRAM  I 

SURFACE  OF  HELICAL  PINION  GENERATED  BY  CONE  CUTTER 

AUTHORS:  FAYDOR  LITVIN 
JIAO  ZHANG 


* 

* 

* 

* 

* 

* 

* 

* 


******************************************************* 


PURPOSE 

THIS  PROGRAM  IS  USED  TO  CALCULATE  THE  SURFACE  OF  A  HELICAL  PINION 
WHICH  IS  GENERATED  BY  CONE  CUTTER 


NOTE 


THIS  PROGRAME  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY 
COMPILE  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,KSIC,MUO 

DIMENSION  X ( 1 00) , Y ( 1 00) . ERROR ( 1 00) , FEEREC (100), UREC (100), 
+  THETAR (100) 


DEFINE  PARAMETERS  USED  BY  PROGRAMS 

(1)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INPUT  AND  OUTPUT 
DEVICES 

IN=5 
LP=6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG=1 

C  (3)  NSOLVE  IS  THE  UPPER  LIMITATION  OF  REPEATATION  FOR  SOLVING 
C  NONLINEAR  EQUTIONS; 

C  EPSI  IS  THE  CLEARANCE  OF  FUNCTION  VALUES  WHEN  THE  FUNCTIONS 

C  IS  CONSIDERED  AS  SOLVED  (ALL  FUNTIONS  HAVE  FORMS  OF  F(X)=0); 

C  DELTA  IS  THE  CLERANCE  OF  VARIABLE  INCREMENT  WHEN  FUNCTION  IS 

C  SOLVED 

C  NC,  EPSI  AND  DELTA  MAY  BE  CHANGED  WHEN  SOLUTIONS  ARE  DIVERGENT 

C  OR  LESS  ACCURATE 

NSOLVE'lOO 
DELTA' l.D- 15 
EPSI'l .D~15 

C  (A)  OTHER  PARAMETERS (DON ' T  CHANGE) 

DR'DATAN ( 1 . DO) / A5 . DO 
C 

C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM  (USE  INCH  AS  UNIT  OF  LENGTH) 

C 

C  (1)  PINION  AND  GEAR:  PN-DIAMETRAL  PITCH;  Nl-PINION  TOOTH  NUMBER; 

C  KSIN-PRESSURE  ANGLE  IN  NORMAL  SECTION  (DEGREE) ; 

C  BETAP-HELIX  ANGLE  (DEGREE) ; 

C  HD=HEIGHT  OF  DEDENDUM  OF  PINION; 


95 


C  HA-HEIRHT  OF  ADDENDUM  OF  PINION 

C  ZCOE-COEFFICIENT  OF  PINION  TOOTH  LENGTH  (THE  LENGTH-  ZCOE/PN) 

PN-10.D0 
Nl-20 

KS IN-20. D0*DR 
BETAP-20 . D0*DR 
HD-l.DO/PN 
HA- 1 .  DO/ PN 
ZCOE-IO. 

C  (2)  TOOL:  ALPHA-HALF  OF  CONE  VERTEX  ANGLE (DEGREE) ; 

C  RC-RADIUS  OF  BOTTOM  CIRCLE  OF  CONE; 

C  MU-TILT  ANGLE  TO  INSTALL  PINION-CUTTING  TOOL 

ALP-80. D0*DR 

MUO-DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU-1 . *MUO 
RC-1.D0 

C  (3)  OUTPUT: 

C  NZ-NO.  OF  CROSS  SECTIONS  WHERE  PINION  PROFILE  IS  SIMULATED; 

C  NU-NO.  OF  MAXIMUM  POINTS  USED  FOR  SIMULATE  PINION  PROFILE 

C  UIN-INCREMENT  OF  PINION  TOOTH  SURFACE  COORDINATE  U 

NZ-21 
NU-61 

UIN-2.5D0*CH/ (NU-1) 

C 

C  DESCRIPTION  OF  OUTPUT  VARIABLES 

C  Zl-DISTANCE  BETWEEN  CROSS  SECTION  CONSIDERED  AND  MIDDLE  CROSS 
C  SECTION 

C  NO-OUTPUT  NO. 

C  U-TOOL  SURFACE  COORDINATE 

C  THETA-TOOL  SURFACE  COORDINATE 

C  FEE-PINION  TOOTH  SURFACE  GENERATION  PARAMETER 

C  Xl-X  COORDINATE  OF  PINION  PROFILE 

C  Yl-Y  COORDINATE  OF  PINION  PROFILE 

C  Rl-RADIUS  OF  PINION  PROFILE 

C  VSH-AVERAGE  DEVIATION  SHIFT  OF  CROSS  SECTION  PROFILE  FROM  PROFILE 
C  OF  GENERAL  PINION  (INVOLUTE  CURVE) 

C  VPE-MAXIMUM  DEVIATION  OF  CROSS  SECTION  PROFILE  FROM  INVOLUTE  CURVE 
C  VSD-STANDARD  DEVIATION  OF  CROSS  SECTION  PROFILE  FROM  INVOLUTE 

C  CURVE 

C 

C  THE  PROGRAM  IS  WRITTEN  BY  JIAO  ZHANG 
SK-DSIN (KSIN) 

CK-DCOS (KSIN) 

SB-DS IN (BETAP) 

CB-DCOS (BETAP) 

SM-DSIN (MU) 

CM-DCOS (MU) 

SA-DS IN (ALPHA) 

CA-DCOS (ALPHA) 

KS IC-DATAN (SK/ CK/ CB) 

CKC-DCOS (KSIC) 

SKC-DSIN (KSIC) 

PT-PN*CB 

RP-N1/2./PT 
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RB-RP*DCOS (KSIC) 

RA-RP+HA 

CL-RC/SA 

CH-HD/CK/CM 

Al-CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 

A2»SK*SM*CB-CM*SB 

A3“CA*SK*CM*CB-SA*CK*CB+CA*SM*SB 

A4=SK*CM*CB+SM*SB 

B1-CA*SK-SA*CK*CM 

B2=CK*SM 

B3-SA*SK+CA*CK*CM 

B4-CK*CM 

C1“CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 

C2=SK*SM*SB+CM*CB 

C3-CA*SK*CM*SB-SA*CK*SB-CA*SM*CB 

C4—SK*CM*SB+SM*CB 

D1=CM*CB+SK*SM*SB 

D2=SA*SM*CB-SB* (CA*CK+SA*SK*CM) 

D3=CA*CK*SB 

D4-CA* (ca*sk-sa*ck*cm) 

D5=SA* (sa*sk+ca*ck*cm) 

D6=CA*CK*SM 
DO  5  1=1, NZ 

Zl— ZCOE/PN/2.+ZCOE/PN*FLOAT(l-l)/FLOAT(NZ-l) 

FEEO=Z 1 *SB/ CB/RP 
CFO=DCOS (FEEO) 

SFO=DSIN (FEEO) 

R1=RP 

ERR=0. 

NP=0 

DO  15  J-l.NU 

IF  (R 1 .GT.RA.AND. J,GT. 1)  GOTO  55 
U=CL-UIN*FLOAT(J-l) 

TEMPI- (Zl-(CL-CH) *C4-U*CA*C3)/U/SA/DSQRT(C1*C1+C2*C2) 
TEMP2=DARSIN(C1/DSQRT(C1*C1+C2*C2)) 

TEMP3-DARS IN (TEMP 1 ) 

THETA-TEMP3-TEMP2 
CT-DCOS (THETA) 

ST-DS IN (THETA) 

XC=U*SA*  (CT*A1+ST*A2)  +U*CA*A3-  (CL-CH)  *A4 
YC=U*SA* (CT*B1-ST*B2) -U*CA*B3+ (CL-CH) *B4 

FEE- (U* (CT*Dl+ST*D2) - (CL-CH) * ( (CT*CA*CA+SA*SA) *D1~ST*D3) ) /RP/ 
+  (CT*D4+D5-ST*D6) 

CF-DCOS (FEE) 

SF-DSIN (FEE) 

X1=XC*CF+YC*SF-RP*FEE*CF+RP*SF 
Yl— XC*SF+YC*CF+RP*FEE*SF+RP*CF 
Rl-DSQRT (X1*X1+Y1*Y1) 

IF  (Rl.LT.RB)  GOTO  15 

NP-NP+1 

X(NP)-X1 

Y  (NP) -Yl 

FEEREC (NP)-FEE 

UR EC (NP)-U 
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JL 


V 


' 


* 


THETAR (NP) -THETA 
XE-X1*CF0+Y1*SF0 
YE— X1*SF0+Y1*CF0 
YS-YE 

FEET-FEE+FEEO 
DO  25  K-l , NSOLVE 

w-  RP*FEET*CKC*DS IN (FEET-KSIC) +RP*DCOS (FEET) -YS 
DWDF— RP*SKC*DCOS (FEET-KSIC) +RP*FEET*CK*DCOS (FEET-KSIC) 

DF— W/DWDF 

IF  (DABS (W) .LT. EPS I. AND. DABS (DF) .LT. DELTA)  GOTO  65 
FEET-FEET+DF 
25  CONTINUE 
65  CONTINUE 

XS— RP*FEET*CKC*DCOS (FEET-KSIC) +RP*DSIN (FEET) 

ERROR  (NP)-XE-XS 

C  WRITE  (LP.110)  K, W,DF, ERROR (NP) ,FEET,FEE+FEEO 
C  110  FORMAT  (IX, ’K= ' , 15 , 5X, ' W= ’ ,D15 . 7, 5X, ' DF= 1 ,D15 . 7,D15 . 7, 2D25 .17) 
ERR-ERR+ERROR (NP) 

15  CONTINUE 
55  VSH— ERR/FLOAT  (NP) 

VPE-O. 

VSD-O. 

DO  35  J-l.NP 
ETEMP-ERROR (j) +VSH 

IF  (DABS (ETEMP) . GT . VPE)  VPE-DABS (ETEMP) 

VSD-VSD+ETEMP**2 
35  CONTINUE 

VSD=DSQRT (VSD/FLOAT  (NP) ) 

WRITE  (LP , 10)  Z1  , VSH, VPE.VSD 

10  FORMAT (1H1, ///' Zl= ' ,F15. 7, 5X, ' VSH- 1 , F15 . 7 , 5X, ’ VPE= ' , F15 . 7 , 5X, 

+  ' VSD' ,D15. 7.5X/2X, 'NO. ’ ,7X, 'U' ,1AX, ’THETA' ,10X, 'FEE' ,12X, 

+  ’XI ' ,13X, ' Y1 ’ ,13X, ' R1 ’ ,13X, 'ERROR') 

DO  A5  J-l.NP 

Rl-DSQRT (X ( J) **2+Y ( J) **2 ) 

ETEMP-ERROR ( J) +VSH 

WRITE  (LP.20)  J,UREC(J) ,THETAR(J) ,FEEREC(J) ,X(J) ,Y(J) ,R1, ETEMP 
20  FORMAT  (2X, 12, 2X, 7F15 . 7) 

45  CONTINUE 
5  CONTINUE 
STOP 
END 

SENTRY 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 
RP-FLOAT (Nl) /2. DO/PN 
CL-RC/DSIN (ALP) 

D-CL*DCOS (ALP) 

Al-CL-HDC/DCOS (RKS) 

RPU-RP+HAC 
DO  5  1-1 ,NL 
Z-FLOAT (I— 1) *ZI 
YY-D-Z/DTAN (ALP) 

WRITE  (LP, 10)  Z 

10  FORMAT  (1H1/1X, 'Zl- ' .F15.7/1X, 'NO’ ,10X, 'Yl ' ,13X,  'XP'  ,13X,  'YP'  , 

+  13X, 'Rl' ,13X, ’FEE’) 


Jt 


_ #...  ... 
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C  CALCULATE  THE  PROPILE  OF  THE  SURFACE  CUT  BY  PLANE  Z-CONST 
KKK-0 

DO  15  J-1 ,N 

Yl-YY*FLOAT(J-l) /FLOAT (N-l) 

F-DSQRT (1 . - (Z/ (D-Yl) /DTAN (ALP) ) **2) 

FEE* ( (D-Yl) *F/DCOS (ALP) -Al* (F*DCOS (ALP) *DCOS (ALP) +DSIN (ALP) *DSIN ( 
+  ALP) ) ) /RP/ (DSIN (ALP) *DCOS (ALP-RKS) -F*DCOS (ALP) *DSIN (ALP-RKS)  ) 
XP (I , J) - (D-Yl) * (DTAN (ALP) *F*DCOS (ALP-RKS+FEE) -DSIN (ALP-RKS+FEE) ) 

+  +Al*DSIN (FEE-RKS) -RP*FEE*DCOS (FEE) +RP*DSIN (FEE) 

YP (I ,  J)  =- (D- Yl)  * (DTAN (ALP) *F*DSIN (ALP-RKS+FEE) +DCOS (ALP-RKS+FEE) ) 
+  +A1 *DCOS (FEE-RKS) +RP*FEE*DS IN (FEE) +RP*DCOS (FEE) 

R1=DSQRT(XP (I , J) **2+YP (I , J) **2) 

IF  (Rl.GT.RPU)  THEN 
JJ-J-1 

XP (I , J) -XP (I , JJ) + (XP (I , J) -XP (I , JJ) ) / (R1-R1TEMP) * (RPU-R1TEMP) 

YP (I , J) =YP (I , J J) + (YP (I , J) -YP (I , J J) ) / (R1-R1TEMP) * (RPU-R1TEMP) 
FEE-FEETEM+ (FEE-FEETEM) / (Rl-RITEMP) * (RPU-R1TEMP) 

Y1*Y1TEM+ (Y1-Y1TEM) / (Rl-RITEMP) * (RPU-R1TEMP) 

Rl*DSQRT (XP (I , J) **2+ YP (I , j) **2) 

KKK=1 
END  IF 

WRITE  (LP , 20)  J,Y1,XP(I,J) , YP (I , J) , Rl , FEE 
20  FORMAT  (1X,I4,4F15.7,F15.7) 

IF  (KKK.GT.O)  GO  TO  30 
FEETEM=FEE 
Y1TEM=Y1 
R1TEMP=R1 
15  CONTINUE 
30  NS (i) “J~l 

IF  (I.NE.l)  GO  TO  55 
NSl-NS(l) 

GO  TO  5 

C  PREPARATION  OF  INTERPLORATION 
55  NS2-NS(I) 

DO  105  L=1 ,NS2 
J=NS2+1-L 

IF  (YP(1,NS1) .GT.YP(I, J))  GO  TO  110 
105  CONTINUE 
110  NS2-J 
NREC=2 
XERS-O.DO 
DO  115  L=1,NS2 
DO  125  J-NREC.NS1 
IF  (YP  (1 ,  J)  .GT.  YP  (I ,  L)  )  GO  TO  120 
125  CONTINUE 
120  NREC-J 
Jl-J-1 

XERROR  (L)  -XP  (I ,  L)  -XP  (1 ,  Jl)  -  (YP  (I ,  L)  -YP  (1 ,  Jl)  )  *  (XP  (1 ,  J)  -XP  (1 ,  Jl)  ) 

+  /(YP(1,J)-YP(1,J1)) 

XERS-XERS+XERROR (L) 

115  CONTINUE 

XERS— XERS/FLOAT  (NS2) 

XPE-O.DO 

SDX-O.DO 
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IF  (NDBUG.GT. 2)  WRITE  (LP,80) 

80  FORMAT  (1H1.13X, 'NO. ' ,4X, 'DIVIATION  VALUE') 

DO  45  L-1.NS2 

XERROR (L) “XERROR (L) +XERS 

IF  (NDBUG.GT. 2)  WRITE  (LP,50)  L, XERROR (L) 

50  FORMAT  (13X,I3,2F15.7) 

IF  (DABS (XERROR (L) ) . GT . XPE)  XPE-DABS (XERROR (L) ) 
SDX=SDX+XERROR (L) **2 
45  CONTINUE 

SDX-DSQRT (SDX/FLOAT (NS2) ) 

WRITE  (LP.40)  XERS ,XPE , SDX 

40  FORMAT  (///IX, ' XSH= ' ,E15.7,5X, 'XPE=' ,E15.7,5X, 'SDX“' ,E15.7) 
IF  (NDBUG.GT. 2)  WRITE  (LP,60)  NS1.NS2 
60  FORMAT  (//5X, ’ NS1- ’ ,I6,6X, 'NS2-’ ,16) 

5  CONTINUE 
STOP 
END 
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PROGRAM  II 

UDDERCUTTING  CONDITION  FOR  HELICAL  PINION 
GENERATED  BY  CONE  CUTTER 
AUTHOR:  FAYDOR  LITVIN 
JIAO  ZHANG 


* 

* 

* 

* 

* 

* 

* 

* 


******************************************************* 


PURPOSE 

THIS  PROGRAM  IS  USED  TO  FIND  THE  UDDERCUTTING  CONDITIONS  FOR  A 
HELICAL  PINION  GENERATED  BY  CONE  CUTTER 


NOTE 


THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY 
COMPILE  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM. 

IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN.MU 


DEFINE  PARAMETERS  USED  BY  PROGRAMS 


(1)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INPUT  AND  OUTPUT 
DEVICES 

IN-5 
LP-6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG=2 

C  (3)  OTHER  PARAMETERS (DON ' T  CHANGE) 

DR-DATAN (1 . DO) /45 . DO 
C 

C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM (USE  INCH  AS  UNIT  OF  LENGTH) 

C 

C  (1)  PINION  AND  GEAR:  PN=DIAMETRAL  PITCH;  Nl-PINION  TOOTH  NUMBER; 

C  KSIN-PRESSURE  ANGLE;  BETAP-HEuIX  ANGEL  OF  PINION; 

C  HD-HEIGHT  OF  DEDENDUM  OF  PINION 

PN-10.D0 
Nl-20 

KS IN-20. D0*DR 

BETAP-30.D0*DR 

HD-1.D0/PN 

C  (2)  TOOL:  ALPHA-HALF  OF  CONE  VERTEX  ANGLE (DEGREE) ; 

C  RC-RADIUS  OF  BOTTOM  CIRCLE  OF  CONE; 

MU-TILT  ANGLE  OF  MOUNTING  TOOL 
ALPHA- 20. D0*DR 
RC-0.35D0 

MUO-DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU-0 . *MU0 

C  (3)  PROBLEM:  NPROB-ID  NO.  OF  PROBLEM  (-1-GIVEN  N1  AND  HD,  FIND  IF 
C  UDDERCUTTING  OCCUR;  0-GIVEN  Nl,  FIND  MAXIMUM  HD  WITHOUT  UDDER- 
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non 


CUTTING;  1-GIVEN  HD,  FIND  MINIMUM  N1  WITHOUT  UDDERCUTTING); 
N-NO.  OF  THETA  VALUES  USED  CALCULATION; 

DELTHE- INCREMENT  OF  THETA (DEGREE) 

NPROB-  1 
N-21 

DELTHE- l.D0*DR 
C 

C  DESCRIPTION  OF  OUTPUT 
C 

C  OUTPUT  IS  A  STATEMENT  BASED  ON  THE  PROBLEM  WITHOUT  ANY  LITERAL 
C  PARAMETER 
C 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 
SK-DSIN (KSIN) 

CK-DCOS (KSIN) 

SB-DSIN (BETAP) 

CB-DCOS (BETAP) 

SM-DSIN (MU) 

CM-DCOS (MU) 

SA-DS IN (ALPHA) 

CA-DCOS (ALPHA) 

PT-PN*CB 

RP-N1/2./PT 

CL-RC/SA 

CH-HD/CK/CM 

UU-CL-CH 

AA=CM*CB+SK*SM*SB 

BB=CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 

CC=CA*CK*SB 

DD1=CA*SK-SA*CK*CM 

DD=DD1*CA 

EE 1 =SA*SK+CA*CK*CM 

EE-EEl*SA 

FF1=CK*SM 

FF=FF1*CA 

II- (CA*SK*CM*SB-SA*CK*SB-CA*SM*CB) *CA/SA 
WRITE  (LP,90) 

90  FORMAT (1H1) 

IF  (NPROB)  5,15,25 
C  CHECK  IF  UNDERCUTTING  OCCURS 
5  WRITE  (LP, 10) 

10  FORMAT  (1H1/3X, 'NO' ,7X, 'THETA' ,12X, 'A' ,13X, * B ’ ,14X, ' C ' , 10X, 

+  'B**2-4.*A*C' ,3X, 'U1/(RC/S IN (ALPHA)) ') 

UMIN-5.D0 
DO  45  1=1, N 
NN- (N+l)/2 

THE-DBLE (FLOAT (I-NN) ) *DELTHE 
ST-DSIN(THE) 

CT-DCOS (THE) 

GG-DD*CT+EE-FF*ST 

WW- (BB*BB+AA*AA) *EE+ (BB*DD-AA*FF) *1 I+ST*ST*ST* (AA*AA*FF-BB*BB*FF 
+  +2 . *AA*BB*DD) -CT*CT*CT* (BB*BB*DD~AA*AA*DD+2 . *AA*BB*FF) 

+  -ST* (2 . *AA*AA*FF+AA*BB*DD-AA*EE*I I) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+  +BB*EE*II) 


XX- ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA"AA*EE*CA*CA) *ST+ (AA*FF*CA 
+  *CA-CC*DD))*(BB*CT+AA*ST+II) 

YY-DD1 *SA*CT-FF 1 *SA*ST-EE 1*CA 

ZZ-CM*CK 

A-YY*WW 

B=UU* (WW*ZZ+XX*YY) +RP*GG*GG*GG 

C-UU*UU*XX*ZZ 

D-B*B~4.*A*C 

IF  (D.GE.O.)  THEN 

Ul« (-B-DSQRT (D) )/2./A 

U2- (-B+DSQRT (D) ) /2 . /A 

U1-U1/CL 

IF  (UMIN.GT.U1)  UMIN-U1 

U2-U2/CL 

THE-THE/DR 

IF  (NDBUG.LT. 1)  WRITE  (LP.lOO)  I,THE,A,B,C,D,U1 
100  FORMAT  (IX, I4.8F15. 7) 

ELSE 

IF  (NDBUG.LT. 1)  WRITE  (LP.100)  I , THE, A, B,C,DD 
END  IF 
45  CONTINUE 

IF  (UMIN.LT. 1 .DO)  WRITE  (LP.110) 

110  FORMAT  (///IX, 'UDDERCUTTING  WILL  OCCUR  FOR  YOUR  DESIGN') 

IF  (UMIN.GE.l.DO)  WRITE  (LP,120) 

120  FORMAT  (///IX, 'UDDERCUTTING  WILL  NOT  OCCUR  FOR  YOUR  DESIGN') 

GO  TO  35 

C  DETERMINE  THE  MAXIMUM  ADDENDUM  HEIGHT  OF  RACK  CUTTER 
15  WRITE  (LP, 20) 

20  FORMAT  (1H1/3X, 'NO' ,7X, 'THETA' ,12X, 'A' ,13X, * B ' ,14X, ' C ’ ,10X, 

+  'b**2-4.*A*C' ,3X, 'ALLOWED  RATIO  OF  HD/(1/PN)') 

UMIN-5.D0 
DO  55  1=1, N 

THE-DBLE (FLOAT (I-NN) ) *DELTHE 
ST-DSIN (THE) 

CT=DCOS (THE) 

GG=DD*CT+EE-FF*ST 

WW= (BB*BB+AA*AA) *EE+ (BB*DD“AA*FF) *II+ST*ST*ST* (AA*AA*FF~BB*BB*FF 
+  +2 .  *AA*BB*DD)  -CT*CT*CT*  (BB,l'BB*DD-AA*AA,'cDD+2 .  *AA*BB*FF) 

+  -ST* (2 . *AA*AA*FF+AA*BB*DD-AA*EE*I I) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+  +BB*EE*II) 

XX= ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA-AA*EE*CA*CA) *ST+ (AA*FF*CA 
+  *CA-CC*DD))*(BB*CT+AA*ST+II) 

YY-DD1 *SA*CT~FF 1 *SA*ST_EE 1 *CA 

ZZ-CM*CK 

A-XX*ZZ 

B-CL* (WW*ZZ+XX*YY) 

C-CL*CL*YY*WW+CL*RP*GG*GG*GG 

D-B*B-4.*A*C 

TEST-DABS (A/B) 

EPS-1. D-16 
THE-THE/DR 
IF  (D.GE.O.)  THEN 
IF  (TEST.GT.EPS)  THEN 
U1 - (-B+DSQRT (D) ) /2 . /A 
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U2= (-B-DSQRT (D) ) /2 . /A 

ELSE 

Ul— C/B 

U2-0. 

END  IF 

U1«(CL-U1)*CK*PN 
U2- (CL-U2) *CK*PN 

IF  (NDBUG.LT. 1)  WRITE  (LP.100)  I,THE,A,B,C,D,U1,U2 
ELSE 

IF  (NDBUG.LT. 1)  WRITE  (LP.lOO)  I,THE,A,B,C,D 
END  IF 

IF  (UMIN.GT.U1)  UMIN-U1 
55  CONTINUE 

WRITE  (LP.200)  UMIN 

200  FORMAT  (///IX, 'TO  AVOID  UDDERCUTTING,  IT  IS  NECESSARY  TO  KEEP  DEDE 
+NDUM  OF  PINION  <=' ,F10.7, '/PN') 

GO  TO  35 

C  DETERMINE  THE  MINIMUM  NO.  OF  TEETH  FOR  UNUNDERCUTTING 
25  WRITE  (LP,30) 

30  FORMAT  (1H1/3X, 'NO' ,7X, 'THETA' ,7X, 'NO.  OF  TEETH') 

RNMAX-0. 

DO  65  1=1, N 

THE-DBLE (FLOAT (I-NN) ) *DELTHE 
ST=DSIN (THE) 

CT-DCOS (THE) 

GG=DD*CT+EE-FF*ST 

WW= (BB*BB+AA*AA) *EE+ (BB*DD~AA*FF) *II+ST*ST*ST* (AA*AA*FF~BB*BB*FF 
+  +2. *AA*BB*DD)-CT*CT*CT* (BB*BB*DD-AA*AA*DD+2. *AA*BB*FF) 

+  -ST* (2 . *AA*AA*FF+AA*BB*DD-AA*EE*I I) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+  +BB*EE*II) 

XX- ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA~AA*EE*CA*CA) *ST+ (AA*FF*CA 
+  *CA-CC*DD))*(BB*CT+AA*ST+II) 

YY-DD1 *SA*CT-FF 1 *SA*ST"EE 1 *CA 
ZZ=CM*CK 

RR— (CL*CL*YY*WW+CL*UU* (WW*ZZ+XX*YY) +UU*UU*XX*ZZ) /CL/GG**3 

RN=2.*RR*PT 

THE-THE/DR 

IF  (RNMAX.LT. RN)  RNMAX-RN 
IF  (NDBUG.LT. 1)  WRITE  (LP.100)  I , THE , RN 
65  CONTINUE 

WRITE  (LP,300)  RNMAX 

300  FORMAT  (///  IX, 'WITHOUT  UDDERCUTTING,  MINIMUM  TOOTH  NO.  OF  PINION 
+IS: ' ,F11.7) 

35  WRITE  (LP.400) 

400  FORMAT(lHl) 

STOP 

END 
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FLOWCHART  FOR  PROGRAM  III 


****************************** ******* ****** ************ 


*  * 

*  PROGRAM  III  * 

*  CONTACT  ELLIPSIS  FOR  HELICAL  PINION  GENERATED  BY  * 

*  CONE  CUTTER  IN  MESHING  WITH  REGULAR  HELICAL  GEAR  * 

*  * 

*  AUTHORS:  FAYDOR  LITVIN  * 

*  JIAO  ZHANG  * 

*  * 

*  * 


******************************************************* 


C  PURPOSE 

c 

C  THIS  PROGRAM  IS  USED  TO  FIND  THE  SHAPE  AND  ORIENTAION  OF  THE  CONTACT 
C  ELLIPSE  WHEN  A  HELICAL  PINION  CROWNED  BY  CONE  CUTTER  IS  IN 

C  MESHING  WITH  A  REGULAR  HELICAL  GEAR 

C 


C  NOTE 


C 

C  THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY  V 
C  COMPILER  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM. 


C 

IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,MUO,MUT,KSIC,KF,KH,KSP,KQP,KSG,KQG 
DIMENSION  CMM(3,4) ,EFD(3) ,EHD(3) ,RND(3) ,RD(4) ,EFF(3) ,EHF(3) , 

+  RNF (3) ,RC (3) ,RF (3) ,Rl (3) 

C 

C  DEFINE  PARAMETERS  USED  BY  PROGRAMS 
C 

C  (1)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INlUT  AND  OUTPUT 
C  DEVICES 

IN-5 
LP-6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG-2 

C  (3)  OTHER  PARAMETERS (DON ' T  CHANGE) 

DR-DATAN ( 1 . DO) / 45 . DO 
C 

C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM  (USE  INCH  AS  UNIT  OF  LENGTH) 

C 

C  (1)  PINION  AND  GEAR:  PN-DIAMETRAL  PITCH;  Nl-PINION  TOOTH  NUMBER; 
C  MPG-TOOTH  NUMBER  RATIO (GEAR  TOOTH  NO./N1); 

C  KSIN-PRESSURE  ANGLE  IN  NORMAL  SECTION; 

C  BETAP-HELIX  ANGEL; 

C  /.^-HEIGHT  OF  DEDENDUM  OF  PINION 

PN-10.D0 
Nl-20 
MPG-2 

KSIN-20.D0*DR 
BETAP- 1 0 . DO*DR 
HD-l./PN 

C  (2)  TOOL:  ALPHA-HALF  OF  CONE  VERTEX  ANGLE (DEGREE) ; 

C  RC-RADIUS  OF  BOTTOM  CIRCLE  OF  CONE; 
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n  o 
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C  MU-TILT  ANGEL  TO  INSTALL  PINION  CUTTING  TOOL 

ALPHA-80. DO*DR 

MUO-DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU-0 . *MUO 
RC-l.DO 

C  (3)  DEFORMATION:  DEL-CONTACT  DEFORMATION  AT  CONTACT  POINT 
DEL-4. D-4 

(4)  OUTPUT:  NU-NUMBER  OF  CONTACT  POINTS  IN  MATING  SURFACES  FOR  US 
TO  CALCULATE  CONTACT  ELLIPSES) 

NU-101 
C 

C  DESCRIPTION  OF  OUTPUT  PARAMETER 
C 

C  Rl-PINION  RADIUS  OF  CONTACT  POINT 

C  ALPHA-THE  ROTATION  ANGLE  BETWEEN  PRINCIPAL  DIRECTION  OF  PINION 
C  TOOTH  SURFACE  AND  AXES  OF  CONTACT  ELLIPSE 

C  A-LENGTH  OF  HALF  SHORT  AXIS  OF  CONTACT  ELLIPSE 
C  B-LENTH  OF  HALF  LONG  AXIS  OF  CONTACT  ELLIPSE (ALONG  DIRECTION  OF 
C  GEAR  TOOTH  LEHGTH) 

C  RNF-UNIT  NORMAL  OF  PINION  TOOTH  SURFACE  AT  CONTACT  POINT 
C  EFF-PRINCIPAL  DIRECTION  OF  PINION  TOOTH  SURFACE  AT  CONTACT  POINT 

C  EHF-PRINCIPAL  DIRECTION  OF  PINION  TOOTH  SURFACE  AT  CONTACT  POINT 

C 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 
SK-DSIN (KSIN) 

CK-DCOS (KSIN) 

SB-DSIN (BETAP) 

CB-DCOS (BETAP) 

SM-DSIN (MU) 

CM-DCOS (MU) 

SA-DSIN (ALPHA) 

CA-DCOS (ALPHA) 

KS IC-DATAN (SK/CK/CB) 

CKC-DCOS (KSIC) 

SKC-DSIN (KSIC) 

PT=PN*CB 
RP-N1/2./PT 
RB=RP*DCOS (KSIC) 

RA-RP+l./PN 

N2-N1*MPG 

RG-RP*MPG 

CL-RC/SA 

CH-HD/CK/CM 

A-CL-CH 

CMM(1 , 1) -CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 
CMM ( 1 , 2) -SA*CK*CB-CA*SK*CM*CB-CA*SM*SB 
CMM (1 , 3) -SK*SM*CB-CM*SB 
CMM (1 , 4) — (SK*CM*CB+SM*SB) 

CMM (2,1) “CA*SK-SA*CK*CM 
CMM (2 , 2) -SA*SK+CA*CK*CM 
CMM (2,3) — CK*SM 
CMM (2, 4) -CK*CM 

CMM(3,1)-CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 
CMM (3 , 2) -SA*CK*SB-CA*SK*CM*SB+CA*SM*CB 


CMM  (3 , 3) -SK*SM*SB+CM*CB 
CMM (3 , 4) — SK*CM*SB+SM*CB 
D1-CM*CB+SK*SM*SB 
D2=SA*SM*CB-SB* (CA*CK+SA*SK*CM) 

D3-CA*CK*SB 

D4-CA* (C‘*SK-SA*CK*CM) 

D5-SA* (SA*SK+CA*CK*CM) 

D6-CA*CK*SM 

UL-CL-2.D0*CH 

UU-CL 

MUT-MUO/DR 

WRITE  (LP, 110)  RP,RA,RB,MUT 

110  FORMAT  (///IX, 'RP= ' , F15 . 7 , 5X, ' RA= 1 , F15 . 7 , 5X, 1 RB= ' , F15. 7 , 5X, 

+  'MUO= 1 , F15 . 7) 

TH-O.DO 

ST-DSIN(TH*DR) 

CT-DCOS (TH*DR) 

EFD(1)=ST 
EFD(2)=0.D0 
EFD(3)— CT 
EHD(l)— CT*SA 
EHD (2) -CA 
EHD(3)— ST*SA 
RND(1)=CT*CA 
RND(2)-SA 
RND(3)-ST*CA 

CALL  MATMUL(CMM,EFD,EFF,3,3) 

CALL  MATMUL (CMM , EHD , EHF , 3 , 3) 

CALL  MATMUL(CMM,RND,RNF,3,3) 

WRITE  (LP.10)  TH,  (EFF  (M)  , M-l , 3)  ,  (EHF  (M)  ,M=1 , 3)  ,  (RNF  (M)  ,M=1 , 3) 

10  FORMAT (1H1.///3X, ' THTEA: ' ,F15.7, '  DEGREE' //IX, 'EFF: ’ .3F15.7/1X, 

+  ’EHF: ' .3F15.7/1X, 'NF  :',3F15.7) 

DO  15  J-l.NU 

U=UU- (UU-UL) /FLOAT (NU-1) *FLOAT ( J- 1) 

KF—CA/SA/U 

KH=0.D0 

RD  ( 1 ) =U*CT*SA 

RD(2) ”-U*CA 

RD(3)=U*ST*SA 

RD(4) “A 

CALL  MATMUL (CMM , RD , RC , 3 , 4) 

FEE” (U* (CT*D1+ST*D2) -A* ( (CT*CA*CA+SA*SA) *D1-ST*D3) ) / (CT*D4+D5~ST* 
+  D6)/RP 
RF (1) -RC (1) -RP*FEE 
RF(2)-RC(2)+RP 
RF(3)-RC(3) 

RPP-DSQRT (RF (1) **2+RF (2)  **2) 

CALL  PROT (RF,R1,FEE) 

TEMP-RC (2) *EFF (1) -RF (1)  *EFF (2) 

B13-EHF(3)-KF*TEMP 
B23— EFF  (3) 

B33—  (RNF  ( 1)  *RF  (1)  +RNF  (2)  *RF  (2)  +KF*TEMP*TEMP) 

SIGMA-0. 5D0*DATAN(2.*B13*B23/(B23*B23-B13*B13- (KF-KH) *B33)) 

COE1- (B23*B23-B13*B13- (KF-KH) *B33) /B33/DCOS (2. *SIGMA) 
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C0E2-KF+KH+ (B13*B13+B23*B23) /B33 
KSP= (C0E2-C0E 1 ) /2 . DO 
KQP= (COE2+COE 1 ) /2 . DO 

IF  (NDBUG.LT. 1)  WRITE  (LP,20)  U,FEE,RPP, (R1 (M) ,M=1,3)  , 

+  B13.B23.B33, SIGMA, KSP.KQP 

20  FORMAT  (//IX,  'U  =  1  ,  F15 . 7 , 5X ,  '  FEE*  '  ,  F15 . 7 ,5X,  '  R1  =’,F15.7/ 

+  IX, 'XP  ,F15.7,5X, ’YP  ,F15.7,5X, 'ZP  =’,F15.7/ 

+  IX, ’ B 1 3=  * ,F15.7,5X, 'B23=' ,F15.7,5X, ’B33-’ ,F15.7/ 

+  IX, ' SIG= ' ,F15.7,5X, 1  KSP= ' , F15 . 7 ,5X, ' KQP= ' , F15 . 7) 

W—  (A-U) 

KSG-O. 

KQG-1 . / (RG*SK* (CKC/CK/CB) **2+W*CM*CK/SK) 

G1=KSP-KQP 

G2=KSG-KQG 

S1=KSP+KQP 

S2=KSG+KQG 

SIGMGP=-SIGMA 

ALPHA1=0. 5D0*DATAN (G2*DSIN (2. *SIGMGP) / (G1~G2*DC0S (2. *SIGMGP) ) ) 
ALPHA2=ALPHAl+SIGMGP 

AA= (S1-S2-DSQRT (Gl*Gl-2. *G1*G2*DC0S (2. *SIGMGP) +G2*G2) ) /4. 

BB= (S1-S2+DSQRT (Gl*Gl-2. *G1*G2*DC0S (2. *SIGMGP) +G2*G2) ) /4. 

AAA* 1 . /DSQRT (DABS (AA)  ) 

BBB* 1 . /DSQRT (DABS (BB) ) 

RATIO=BBB/AAA 

ALPHAl=ALPHAl/DR 

ALPHA2=ALPHA2/DR 

C  WRITE  (LP.130)  G1 , G2 , S 1 , S2 , ALPHAl 
C  130  FORMAT  (1X.5F15.7) 

WRITE  (LP.30)  RPP.ALPHA2, AAA, BBB, RATIO 
30  FORMAT  (IX, 'R1  - ' , FIS . 7 , 5X, ' ALP* ' , F15 . 7 , 5X, ' A  »',F15.7,5X, 

+  'B  =',F15.7,5X,’B/A=',F15.7) 

IF  (RPP.GT.RA)  GO  TO  5 
15  CONTINUE 
5  STOP 
END 
C 
C 

SUBROUTINE  MATMUL (CMM , A, B , N , M) 

C  THIS  SUBROUTINE  IS  USED  TO  MULTIPLY  THE  MATRIX  CMM (N*M)  BY  THE  MATRIX 
C  A(M*1).  THE  RESULT  IS  STORED  IN  THE  MATRIX  B(N*1) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  CMM (3 , M)  ,  A (M)  ,  B  (N) 

DO  5  1=1,  N 
5  B(I)=0. 

DO  15  1=1, N 
DO  15  J=1,M 

15  B(I)=B(I)+CMM(I,  J)*A(J) 

RETURN 

END 

C 

C 

SUBROUTINE  PROT (A, B , FEE) 

C  THIS  SUBROUTINE  IS  USED  TO  ROTATE  COORDINATE  SYSTEM  PLANARLY  IN  XOY 
C  THROUGH  ANGLE  FEE 
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DOUBLE  PRECISION  A (3) , B (3) , FEE 
B (1) -A(l) *DCOS (FEE) +A(2) *DSIN (FEE) 
B (2) =A (2) *DCOS (FEE) -A ( 1 ) *DS IN (FEE) 


B(3)-A(3) 

RETURN 

END 
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******************************************************* 


THIS  PROGRAM  IS  USED  TO  CALCULATE  THE  TRANSMISSION  ERRORS  OF  A 
PINION  GENERATED  BY  THE  CUTTER  WITH  CONE  OR  REVOLUTE  SURFACE 
IN  MESHING  WITH  A  MISALIGNED  REGULAR  HELICAL  GEAR 


THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY  V 
COMPILER  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM. 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN.MU.MUO.KSIC 
DIMENSION  Z  (99)  .ANGLE  (99)  ,  ERROR  (99)  ,  ERR  (99) 

COMMON  /BLOCK1/  X(1I) , Y(10) ,A(10, 10) , Y1 (10) , IPVT(IO) .WORK(lO) , 

+  EPS I , DELTA, NC , NE , NDIM 

COMMON  /BLOCK2/  S(4,4) ,CMM(4,4) ,CMCD(4, 4) ,C, R.RP.RG.Al .HD.RL.KSIC, 
+  ALPHA . CK , SK , CB , SB , CM , SM , CA , SA , CKC , SKC , NTOOL 


DEFINE  PARAMETERS  USED  BY  PROGRAMS 


(1)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INPUT  AND  OUTPUT 
DEVICES 

IN-5 
LP-6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG-1 

C  (3)  NC  IS  THE  UPPER  LIMITATION  OF  REPEATATION  FOR  SOLVING  NONLINEAR 

C  EQUTIONS; 

C  EPSI  IS  THE  CLEARANCE  OF  FUNCTION  VALUES  WHEN  THE  FUNCTIONS 

C  IS  CONSIDERED  AS  SOLVED  (ALL  FUNTIONS  HAVE  FORMS  OF  F(X)-O); 

C  DELTA  IS  THE  RELATIVE  DIFFERENCE  FOR  TAKING  DERIVATIVES 

C  NC,  EPSI  AND  DELTA  MAY  BE  CHANGED  WHEN  SOLUTIONS  ARE  DIVERGENT 

C  OR  LESS  ACCURATE 

NC-100 
DELTA-l.D-3 
EPSI-l.D-12 

C  (4)  OTHER  PARAMETERS (DON ' T  CHANGE) 

NDIM- 10 
NE-5 
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DR=DATAN (1 . DO) /45 . DO 

C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM (USE  INCH  AS  UNIT  OF  LENGTH) 

C  (1)  PINION  AND  GEAR:  PN-DIAMETRAL  PITCH;  NP=PINION  TOOTH  NUMBER; 

C  RMPG=TOOTH  NUMBER  RATIO (GEAR  TOOTH  NO./NP); 

C  KSIN“PRESSURE  ANGLE  IN  NORMAL  SECTION; 

C  BETAP=HELIX  ANGLE  OF  PINION  AND  GEAR; 

C  HD=HEIGHT  OF  DEDENDUM  OF  PINION; 

C  COE=COEFF .  OF  CENTRAL  DISTANCE (USUALLY  COE=l.) 

PN=10.D0 

NP=20 

RMPG=2.D0 

KSIN=20.D0*DR 

BETAP=15.DO*DR 

HD=1 .DO/PN 

COE=l . OOODO 

C  (2)  TOOL:  ALPHA=HALF  OF  CONE  VERTEX  ANGLE (DEGREE) ; 

C  RC=RADIUS  OF  BOTTOM  CIRCLE  OF  CONE 

C  R=RADIUS  OF  ARC 

C  MU=TILT  ANGEL  OF  MOUNTING  TOOL 

NTOOL=TOOL  ID  NO.  (l=CONE  SURFACE,  2=REVOLUTE  SURFACE) 

NTOOL=2 

ALPHA=20.D0*DR 
RC=10.3527D0 
R=3.0D1 

MUO=DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU=  0 . *MUO 

(3)  MISALIGNMENT:  NMIS=ID  NO. (l=CROSSING  AXES,  2= INTERSECTING  AXES); 
NG=NO.  OF  MISALIGNED  ANGLES  TO  BE  SIMULATED  (FROM  -(NG-l)/2  TO 
(NG-l)/2  TIMES  GAMMA I  WITH  ODD  NG) ; 

GAMMAI= INCREMENT  OF  MISALIGNED  ANGLE (MINUTE) ; 

NMIS=2 
NG=2 

GAMMAI=5 .  DO 

C  (4)  OUTPUT:  FEE 1= INCREMENT  OF  ROTATION  ANGLE  OF  PINION (DEGREEE) 
FEEI=1 . ODO*DR 
C 

C  DESCRIPTION  OF  OUTPUT  PARAMERTERS 
C 

C  FEEl=ROTATION  ANGLE  OF  PINION 

C  FEE2=ROTATION  ANGLE  OF  GEAR 

C  RP'RADIUS  OF  PINION  CONTACT  POINT 

C  RG=RADIUS  OF  GEAR  CONTACT  POINT 

C 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 

C 

C 

SK=DSIN (KSIN) 

CK=DCOS (KSIN) 

SB-DS IN (BETAP) 

CB=DCOS (BETAP) 

SM=DSIN(MU) 

CM-DCOS (MU) 

SA»DS IN (ALPHA) 

CA-DCOS (ALPHA) 
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KSIC-DATAN (SK/CK/CB) 

ckodcos  (ksic) 

SKC=DSIN (KSIC) 

PT“PN*CB 
RP-NP/2./PT 
RB=RP*DCOS (KSIC) 

RA-RP+l./PN 
NG=NP*RMPG 
RG-RP*RMPG 
RL“RC*CA/SA 
C  CH=HD/CK/CM 

Al=RL/CA-HD/CK/CM 
C=  (RP+RG) *COE 

NCOEF=360 . D0*DR/FEEI /FLOAT (Nl) +0 . 3 

N“NC0EF*2+1 

CALL  INTMAT(CMCD,4,4) 

AA=RL*SA*SA/CA-HD/CK/CM 

CMCD ( 1 , 1)=CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 

CMCD (1 , 2) =SA*CK*CB-CA*SK*CM*CB-CA*SM*SB 

CMCD  (1 , 3)  =‘SK*SM*CB-CM*SB 

CMCD (1,4) — RL*SA*CK*CB-AA* (SK*CM*CB+SM*SB) 

CMCD (2, 1) =CA*SK-SA*CK*CM 

CMCD (2 , 2) =SA*SK+CA*CK*CM 

CMCD (2,3) =-CK*SM 

CMCD (2,4) =-RL*SA*SK+AA*CK*CM 

CMCD (3 , 1) =CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 

CMCD (3 , 2) =SA*CK*SB-CA*SK*CM*SB+CA*SM*CB 

CMCD (3 , 3) =SK*SM*SB+CM*CB 

CMCD (3 , 4) —RL*SA*CK*SB+AA* (-SK*CM*SB+SM*CB) 
CALL  INTMAT (S, 4, 4) 

NGG= (NG-l)/2 
DO  505  LL=1 ,NG 

GAMMA-GAMMAI*FLOAT (LL-NGG) /60 . D0*DR 
CG=DC0S (GAMMA/60 . *DR) 

SG=DS IN (GAMMA/60 . *DR) 

IF  (NMIS.EQ. 1)  THEN 
WRITE  (LP,500)  COE, GAMMA 

500  FORMAT  (IHl , ///, IX, 'C= ' , F4. 2 , (RP+RG) 

+  F5 . 1 ,  ’  (M)  ' ) 

S(1,1)-CG 
S(1,3)=-SG 
S(3, 1)*=SG 
S(3,3)*=CG 
ELSE 

WRITE  (LP,501)  COE, GAMMA 

501  FORMAT  (IHl , /// , IX, ’ C- 1 , F4. 2 , ’ * (RP+RG) 

+  F5.1,’(M)') 

S (2, 2) “CG 
S(2,3)=-SG 
S(3,2)-SG 
S(3,3)-CG 
END  IF 

DO  205  L-1,2 
DO  5  I-l.NE 


CROSSING  ANGLE= ’ , 


INTERSECTING  ANGLE 
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5  X(l)-0.D0 

IF  (NT00L.EQ.1)  X(5)-RL/CA~HD/CK/CM 
C  WRITE  (6,1100)  (X( JK) , JK=1 , 5) ,RL, CA,HD,CK,CM 
C1100  FORMAT  (IX, * ###### ' , 5F15.7/7X.5F15. 7) 

DO  15  I-l.N 

X(7) =FEEI*FLOAT (I- (N+l) /2) 

IF  (L.EQ.l)  X(7)-0.D0 
C  X (5) -DARSIN (RP*X (7) /SK/R) 

X (2) -X (7) /RMPG 
CALL  NONLIN 
X(8)=X(2)+X(3) 

IF  (L.EQ.l)  THEN 
XIN=X(8) 

WRITE  (LP.10) 

10  FORMAT  (////8X, 'FEEl(D) 1 ,8X, ’FEE2(D) ' ,8X, 'K-ERROR(S) ' ,5X, 

+  ' RP' , 13X, 'RG' , F15.7/) 

GO  TO  205 
END  IF 

X(8)=X(8)-XIN 
X(7)“X(7)  /DR 
X(8)=X(8)/DR 

X (9) - (X (8) -X (7) /RMPG) *3600. DO 
Z  (I)  BX  (8) 

ERR  (I)  =X  (9) 

WRITE  (LP , 20)  (X(J) ,J=7,11) 

20  FORMAT  (1X,5F15.7,F15.7) 

WRITE  (LP,30)  (X(J) , J=l,6) 

30  FORMAT  (1X,'#####',6F15.7) 

WRITE  (LP.30)  XP,YP,ZP,XG,YG 
15  CONTINUE 
NT=N-NCOEF 

FII-FEEI/DR/2.D0*FL0AT (ncoef) 

WRITE  (LP.80) 

80  FORMAT  (//,'  FIND  THE  WORKING  RANGE  FOR  ONE  TOOTH: ' , F15 . 7/) 
DO  55  I-l.NT 

X(7)=FEEI*FLOAT(l- (N+l) /2) /DR/2. DO 
X(8)=X(7)+FII 
KK-I+NCOEF 
ANGLE (i) “Z (KK)-Z (I) 

ERROR (I) - (ANGLE (I) -FI I) *3600 . DO 
WRITE  (LP,60)  X(7) ,X (8) .ANGLE (I) .ERROR (I) 

60  FORMAT  (IX, ' (' ,F7.2, ' - ' , F7. 2, ' ) : ' , F15. 7, F15 . 7) 

55  CONTINUE 
DO  95  1*1 , NT 
ATEMP2-ERROR(I) 

IF  (I.NE.l)  THEN 

IF  (ATEMP 1 *ATEMP2 . LE . 0 . DO)  GOTO  105 
END  IF 

ATEMP 1 -ATEMP 2 
95  CONTINUE 

WRITE  (LP, 160) 

160  FORMAT  (//IX, 'MESHING  IS  DISCONTINUOUS') 

105  IF  (DABS (ATEMP 1) .LT. DABS (ATEMP2))  I-I-l 
EMAX-0. 
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EMIN-O. 

DO  135  J-l.NCOEF 
KS-I+J-1 
ET-ERR (KS) 

IF  (ET.LT.EMIN)  EMIN-ET 
IF  (ET.GT.EMAX)  EMAX-ET 
135  CONTINUE 

ET-EMAX-EMIN 

KK-I+NCOEF 

WRITE  (LP, 170)  Z (I) ,Z (KK) ,ET 

170  FORMAT  (//IX, 'WORKING  RANGE  FOR  ONE  TOOTH:  ',F7.2,' - \F7.2/ 

+  IX, 'THE  MAXIMUM  KINEMATIC  ERROR:  ',F15.7,'  (S)',I2) 

205  CONTINUE 
505  CONTINUE 
STOP 
END 

SUBROUTINE  FUNC 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,MUO,KSIC 

COMMON  /BLOCK1/  X (11) , Y (10) , A(10, 10) , Yl  (10) , IPVT (10) , WORK (10) , 

+  EPSI .DELTA, NC,NE,NDIM 

COMMON  /BLOCK2/  S (4, 4) , CMM (4, 4) ,CMCD(4, 4) ,C,R,RP,RG, Al ,HD,RL,KSIC, 
+  ALPHA , CK , SK , CB , SB , CM , SM , CA , SA , CKC , SKC , NTOOL 

DIMENSION  RD(4) ,RND(4) ,RC(4) ,RNC(4) ,RF(4) 

CFEE=DCOS  (X  (3) ) 

SFEE-DSIN (X (3) ) 

CT-DCOS (X (4) ) 

ST-DSIN(X(4)) 

CF-DCOS  (X(7)  ) 

SF-DSIN  (X(7)  ) 

IF  (NTOOL. EQ.l)  THEN 
RD(1)«X(5)  *CT*SA 
RD (2) -RL~X (5) *CA 
RD (3) -X (5) *ST*SA 
RND(1)-CT*CA 
RND(2) -SA 
RND(3)-ST*CA 
ELSE 

CAL-DCOS (ALPHA+X (5) ) 

SAL-DSIN (ALPHA+X (5) ) 

SL2-DSIN(X(5)/2.) 

CAL2-DCOS (ALPHA+X (5) /2 . ) 

SAL2-DSIN (ALPHA+X (5) /2.) 

RD(1)-(A1*SA-2.*R*SL2*SAL2)*CT 

RD (2) -HD/CK/CM*CA+2 . *R*SL2*CAL2 

RD (3) - (A1 *SA-2 . *R*SL2*SAL2) *ST 

RND ( 1 ) -CAL*CT 

RND(2) -SAL 

RND (3) -CAL*ST 

END  IF 

RD(4) -1 . 

CALL  MATMUL (CMCD , RD , RC , 4 , 4) 
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CALL  MATMUL (CMCD , RND , RNC , 3 , 3) 

X (6)  - (RC (1) -RC (2) *RNC (1) /RNC (2) ) /RP 
RF(1) “RC  (1) -RP*X (6) 

RF  (2)  =RC(2)+RP 
RF(3)“RC(3) 

RF (4) “1 . 

CFPF1=DC0S (X (6) +X  (7) ) 

SFPFl-DSIN (X(6)+X(7)) 

XPF=RF (1) *CFPF1+RF (2) *SFPF1 
YPF=RF (2) *CFPF1-RF (1) *SFPF1 
ZPF=RF  (3) 

XNPF=RNC (1) *CFPFl+RNC (2) *SFPFl 
YNPF-RNC (2) *CFPF1-RNC (1) *SFPFl 
ZNPF=RNC (3) 

CF2FG-DC0S  (X(3)) 

SF2FG=DSIN (X(3) ) 

A2=-X(1)*SB+RG*X(2) 

RFG1“CK*CK/CKC*DC0S (X (3) +KSIC) *A2+RG*SF2FG 

RFG2=CK*CK/CKC*DS IN (X (3) +KSIC) *A2~RG*CF2FG 

RFG3=X (1) * (CB+SK*SK*SB*SB/CB) -RG*X (2) *SK*SK*SB/CB 

RNFG1=CK*CB*CF2FG-SK*SF2FG 

RNFG2=CK*CB*SF2FG+SK*CF2FG 

RNFG3“CK*SB 

XGF=RFG1*S (1 , 1) +RFG2*S (1 ,2) +RFG3*S  (1 , 3) 

YGF=RFG1*S (2 , 1) +RFG2*S (2,2) +RFG3*S (2 , 3) 

ZGF=RFG1*S (3,1) +RFG2*S (3 , 2) +RFG3*S (3,3) 

XNGF=RNFG1*S (1 , 1) +RNFG2*S (1 , 2) +RNFG3*S (1 , 3) 

YNGF=RNFG1 *S (2 , 1 ) +RNFG2*S (2 , 2)  +RNFG3*S (2 , 3) 

ZNGF=RNFG1*S (3 , 1) +RNFG2*S (3 , 2) +RNFG3*S  (3 , 3) 

WRITE  (6,100)  X ( 1 )  ,X(2)  ,X(3)  ,X(4)  ,X(5) 

WRITE  (6,100)  XPF , YPF , ZPF , XGF , YGF , ZGF 
100  FORMAT  (IX, ,8F15.7) 

WRITE  (6,100)  XNPF , YNPF , ZNPF , XNGF , YNGF , ZNGF 

Y(1)=XPF-XGF 

Y(2)=YPF-YGF-C 

Y (3) -ZPF-ZGF 

Y (A) =YNPF_ YNGF 

Y (5) “ZNPF-ZNGF 

X (10) =DSQRT (XPF*XPF+YPF*YPF) 

X (1 1) =DSQRT (XGF*XGF+YGF*YGF) 

WRITE  (6,20)  (Y (II) ,11=1,5) 

20  FORMAT  (IX, ' $$$$ ' ,6F15. 7) 

RETURN 

END 

SUBROUTINE  INTMAT  (A,N,M) 

THIS  SUBROUTINE  IS  USED  TO  INITIATE  THE  MATRIX,  WITH  UNIT  DIAGONAL 
ELEMENTS  AND  NULL  OTHER  ELEMENTS 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (A, 4) 

DO  5  I«1,N 
DO  5  J«1,M 
A(I, J)-0. 

IF  (I.EQ.J)  A(I,J)-1. 
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5  CONTINUE 
RETURN 
END 

SUBROUTINE  MATMUL(CMM,A,B,N,M) 

THIS  SUBROUTINE  IS  USED  TO  MULTIPLY  THE  MATRIX  CMCD(N*M)  BY  THE  MATRIX 
A(M*1).  THE  RESULT  IS  STORED  IN  THE  MATRIX  B(N*1 ) 

IMPLICIT  REAL *8 (A~H, 0-Z) 

DIMENSION  CMM (4,4) , A (4) , B (4) 

DO  5  1=1, N 
5  B(I)=0. 

DO  15  1=1, N 
DO  15  J=1,M 

15  B  (I)  “B  (I)  +CMM(l ,  J)  *A(J) 

RETURN 

END 


SUBROUTINE  NONLIN 


IMPLICIT  REAL *8 (A-H, 0-2) 

COMMON  /BL0CK1/  X (1 1) , Y (10) , A(10, 10) , Y1 (10) , IPVT(IO) , WORK(IO) , 
+  EPSI, DELTA, NC,NE,NDIM 

DO  5  1=1, NC 
CALL  FUNC 

WRITE  (6,10)  I, (X(J) ,Y (J) , J=l,5) 

10  FORMAT (IX, '***', 15/5 (IX, 2D15 . 7/) ) 

DO  15  J=1 ,NE 

IF  (DABS (Y (J) ) .GT. EPSI)  GO  TO  25 
15  CONTINUE 
GO  TO  105 
25  DO  35  J=1 ,NE 
35  Yl (J) =Y ( J) 

DO  45  J=1,NE 
DIFF=DABS (X ( J) ) *DELTA 
IF  (X(J) .EQ.O.DO)  DIFF-DELTA 
XMAM=X(j) 

X(J)  =X(J)~DIFF 
CALL  FUNC 
X(J) =XMAM 
DO  55  K*1 ,NE 

55  A(K,  J)  =  (Y1(K)-Y(K))/DIFF 
45  CONTINUE 
DO  65  J=1,NE 
65  Y(J)=-Y1(J) 

CALL  DECOMP  (NDIM,NE,A,COND, IPVT.WORK) 

CALL  SOLVE  (NDIM,NE, A, Y, IPVT) 

DO  75  J=1,NE 
75  X(J)=X(J)+Y(J) 

5  CONTINUE 
105  RETURN 
END 
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SUBROUTINE  DECOMP  (NDIM.N, A, COND, IPVT, WORK) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(NDIM,N) , WORK(N) , IPVT (N) 

C 

C  DECOMPOSES  AREAL  MATRIX  BY  GAUSSIAN  ELIMINATION, 

C  AND  ESTIMATES  THE  CONDITION  OF  THE  MATRIX. 

C 

C  -COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
C  M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL.  1977) 

C 

C  USE  SUBROUTINE  SOLVE  TO  COMPUTE  SOLUTIONS  TO  LINEAR  SYSTEM. 

C 

C  INPUT.. 

C 

C  NDIM  *  DECLARED  ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  A 

C  N  ORDER  OF  THE  MATRIX 

C  A  MATRIX  TO  BE  TRIANGULARIZED 

C 

C  OUTPUT.. 

C 

C  A  CONTAINS  AN  UPPER  TRIANGULAR  MATRIX  U  AND  A  PREMUTED 
C  VERSION  OF  A  LOWER  TRIANGULAR  MATRIX  I-L  SO  THAT 

C  (PERMUTATION  MATRIX) *A=L*U 

C 

C  COND  -  AN  ESTIMATE  OF  THE  CONDITION  OF  A. 

C  FOR  THE  LINEAR  SYSTEM  A*X  =  B  ,  CHANGES  IN  A  AND  B 

C  MAY  CAUSE  CHANGES  COND  TIMES  AS  LARGE  IN  X. 

C  IF  COND+1.0  .EQ.  COND  ,  A  IS  SINGULAR  TO  WORKING 

C  PRECISION.  COND  IS  SET  TO  1.0D+32  IF  EXACT 

C  SINGULARITY  IS  DETECTED. 

C 

C  IPVT  =  THE  PIVOT  VECTOR 

C  IPVT (K)  -  THE  INDEX  OF  THE  K-TH  PIVOT  ROW 

C  IPVT (N)  -  (-1)** (NUMBER  OF  INTERCHANGES) 

C 

C  WORK  SPACE..  THE  VECTOR  WORK  MUST  BE  DECLARED  AND  INCLUDED 
C  IN  THE  CALL.  ITS  INPUT  CONTENTS  ARE  IGNORED. 

C  ITS  OUTPUT  CONTENTS  ARE  USUALLY  UNIMPORTANT. 

C 

C  THE  DETERMINANT  OF  A  CAN  BE  OBTAINED  ON  OUTPUT  BY 
C  DET(A)  =  IPVT (N)  *  A(1 , 1)  *  A (2, 2)  *  ...  *  A(N,N)  . 

C 

IPVT (N) “1 

IF  (N.EQ.l)  GO  TO  150 
NM1-N-1 

C  COMPUTE  THE-  1-NORM  OF  A  . 

ANORM-O.DO 
DO  20  J-1,N 
T-O.DO 
DO  10  I-1,N 
10  T“T+DABS (A(I , J) ) 
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IF  (T.GT.ANORM)  ANORM=T 
20  CONTINUE 


DO  GAUSSIAN  ELIMINATION  WITH  PARTIAL 
PIVOTING. 


DO  70  K=1 ,NM1 
KP1-K+1 


FIND  THE  PIVOT. 


M-K 

DO  30  I-KPl.N 

IF  (DABS (A ( I, K)) .GT.DABS(A(M,K)))  M=I 
30  CONTINUE 
IPVT (K) =M 

IF  (M.NE.K)  IPVT (N)— IPVT (N) 

T-A(M,K) 

A  (M ,  K)  =A  (K ,  K) 

A(K,K) “T 

SKIP  THE  ELIMINATION  STEP  IF  PIVOT  IS  ZERO. 
IF  (T.EQ.O.DO)  GO  TO  70 

COMPUTE  THE  MULTIPLIERS. 

DO  40  I-KPl.N 
40  A(I,K)— A(I,K)/T 

INTERCHANGE  AND  ELIMINATE  BY  COLUMNS. 

DO  60  J-KPl.N 
T=A(M,J) 

A(M,  J)  =A(K,  J) 

A(K, J)“T 

IF  (T.EQ.O.DO)  GO  TO  60 
DO  50  I“KP1 , N 

50  A(I,J)»A(I,J)+A(I,K)*T 

60  CONTINUE 
70  CONTINUE 
C 

C  COND  =  (1-NORM  OF  A) * (AN  ESTIMATE  OF  THE  1-NORM  OF  A~ INVERSE) 

C  THE  ESTIMATE  IS  OBTAINED  BY  ONE  STEP  OF  INVERSE  ITERATION  FOR  THE 
C  SMALL  SINGULAR  VECTOR.  THIS  INVOLVES  SOLVING  TWO  SYSTEMS 
C  OF  EQUATIONS,  (A- TRANSPOSE) *Y  =  E  AND  A *Z  =  Y  WHERE  E 
C  IS  A  VECTOR  OF  +1  OR  -1  COMPONENTS  CHOSEN  TO  CAUSS  GROWTH  IN  Y. 

C  ESTIMATE  «  (1-NORM  OF  Z) / (1-NORM  OF  Y) 

C 

c  SOLVE  (A- TRANSPOSE) *Y  =  E  . 

DO  100  K*1,N 
T-O.DO 

IF  (K.EQ.l)  GO  TO  90 
KM1-K-1 
DO  80  I-l.KMl 
80  T-T+A(I,K)*W0RK(I) 

90  EK= 1 . DO 

IF  (T.LT.O.DO)  EK— 1.D0 
IF  (A(K,K) .EQ.O.DO)  GO  TO  160 
100  WORK 00  —  (EK+T)/A(K,K) 

DO  120  KB-l.NMl 
K-N-RB 
T-O.DO 
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KP1=K+1 

DO  110  I=KP1 ,N 
T=T+A(l,K)*WORK(K) 

WORK (K) =T 
M=IPVT(K) 

IF  (M.EQ.K)  GO  TO  120 
T=WORK (M) 

WORK (M) =WORK(K) 

WORK(K)=T 
120  CONTINUE 
C 

YNORM=O.DO 
DO  130  1=1, N 

130  YNORM=YNORM+DABS (WORK (I) ) 

SOLVE  A *Z  =  Y 

CALL  SOLVE  (NDIM, N, A, WORK, IPVT) 

ZNORM=0.D0 
DO  140  1=1, N 

140  ZNORM=ZNORM+DABS(WORK(I)) 

ESTIMATE  THE  CONDITION. 

COND=ANORM*ZNORM/ YN0RM 
IF  (COND.LT. 1.D0)  COND=l.D0 
RETURN 

C  1-BY-l  CASE.. 

150  COND=l -DO 

IF  (a(1, 1) .NE.O.DO)  RETURN 
C 

C  EXACT  SINGULARITY 

160  COND=l . 0D32 
RETURN 
END 

SUBROUTINE  SOLVE  (NDIM, N, A, B , IPVT) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (NDIM, N) , B (N) , IPVT (N) 

C 

C  SOLVES  A  LINEAR  SYSTEM,  A*X  =  B 

C  DO  NOT  SOLVE  THE  SYSTEM  IF  DECOMP  HAS  DETECTED  SINGULARITY. 

C 

C  -COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
C  M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL,  1977) 

C 

C  INPUT.. 

C 

C  NDIM  =  DECLARED  ROW  DIMENSION  OF  ARRAY  CONTAINING  A 
C  N  ORDER  OF  MATRIX 

C  A  TRIANGULAR IZED  MATRIX  OBTAINED  FROM  SUBROUTINE  DECOMP 

C  B  RIGHT  HAND  SIDE  VECTOR 

C  IPVT  =  PIVOT  VECTOR  OBTAINED  FROM  DECOMP 

C 

C  OUTPUT.. 
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B  -  SOLUTION  VECTOR,  X 


DO  THE  FORWARD  ELIMINATION. 

IF  (N.EQ.l)  GO  TO  50 
NM1-N-1 
DO  20  K=1 ,NM1 
KP1-K+1 
M=IPVT (K) 

T=B (M) 

B(M)-B(K) 

B (K) =T 

DO  10  I=KP1 ,N 
10  B(l)=B(l)+A(l,K)*T 

20  CONTINUE 

C  NOW  DO  THE  BACK  SUBSTITUTION. 

DO  AO  KB=1 ,NM1 
KM1=N-KB 
K=KMl+l 

B(K)=B(K)  /A(K,K) 

T— B  (K) 

DO  30  1=1, KM1 
30  B(I)=B(I)+A(I,K)*T 

40  CONTINUE 
50  B(1)=B(1)/A(1,1) 

RETURN 

END 

C 

C 


l 
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r  **xk 

w*  »  • 

C...  * 

c...  * 
c...  * 
c.. .  * 

C...  *  M 

c...  * 
c...  * 
c...  * 
c...  * 
c...  * 
c...  * 

C . . •  **** 

c 

C  PURPOSE 

C 

C 

C  THIS 
C  CRO' 

C  MOU 

C  AND 

C 

C  NOTE 


★linfcifc********************************^**********^!:*  ****** 


PROGRAM  V  * 

TRANSMISSION  ERRORS  OF  CROWNED  HELICAL  PINION  IN  * 
MESHING  WITH  REGULAR  HELICAL  GEAR  WITH  PREDESIGNED  * 

TRANSMISSION  ERRORS  AND  CONTACT  PATH  * 

* 

AUTHORS:  FAYDOR  LITVIN  * 

JIAO  ZHANG  * 


k  kkk  ******  kkkkkkkkkkkkkkk  kkkkkkkkkkkkkkkkkkkkkkk  kkkkkkk 


THIS  PROGRAM  IS  USED  TO  CALCULATE  THE  TRANSMISSION  ERRORS  OF  A 
CROWNED  HELICAL  PINION  AND  A  HELICAL  GEAR  WHEN  THEIR  AXES  ARE 
MOUNTED  WITH  SOME  ERRORS.  THE  MAGNITUDE  OF  TRANSMISSION  ERRORS 
AND  CONTACT  PATH  ARE  PRE-DESIGNED. 


C 

C  THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY  V 
C  COMPILER  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM. 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,KK,KSIC 

DIMENSION  Z  (99)  .ANGLE (99)  , ERROR  (99)  , ERR  (.99)  ,W(99)  , RPR (99)  ,RGR(99)  , 
+  ZP  (99) ,ZG(99) 

COMMON  /BLOCK1/  X(13) , Y (10) , A(10, 10) , Y1 (10) , IPVT (10) , WORK(IO) , 

+  EPSI, DELTA, NC.NE.NDIM 

COMMON  /BLOCK2/  S(4,4) ,AT(5) ,C,RP,RG,CK,SK,CB,SB,CKC,KSIC,TB, 

+  COEG 1 , COEG2 , RMGP , AA , DR , KK , R 


DEFINE  PARAMETERS  USED  BY  PROGRAMS 


(1)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INPUT  AND  OUTPUT 
DEVICES 

IN=5 
LP-6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG-1 

C  (3)  NC  IS  THE  UPPER  LIMITATION  OF  REPEATATION  FOR  SOLVING  NONLINEAR 
C  EQUTIONS ; 

C  EPSI  IS  THE  CLEARANCE  OF  FUNCTION  VALUES  WHEN  THE  FUNCTIONS 

C  IS  CONSIDERED  AS  SOLVED  (ALL  FUNTIONS  HAVE  FORMS  OF  F(X)=0); 

C  DELTA  IS  THE  RELATIVE  DIFFERENCE  FOR  TAKING  DERIVATIVES 

C  NC,  EPSI  AND  DELTA  MAY  BE  CHANGED  WHEN  SOLUTIONS  ARE  DIVERGENT 

C  OR  LESS  ACCURATE 

NC-100 
DELTA- l.D-3 
EPSI-l.D-10 
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C  (4)  OTHER  PARAMETERS (DON'T  CHANGE) 

NDIM=10 

NE-5 

DR-DATAN ( 1 . DO) /45 . DO 

C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM (USE  INCH  AS  UNIT  OF  LENGTH) 

C  (1)  PINION  AND  GEAR:  PN=DIAMETRAL  PITCH;  NP-PINION  TOOTH  NUMBER; 

C  RMPG-TOOTH  NUMBER  RATIO (GEAR  TOOTH  NO./NP); 

C  KS IN- PRESSURE  ANGLE  IN  NORMAL  SECTION; 

C  BETAP-HELIX  ANGLE  OF  PINION  AND  GEAR; 

C  COE-COEFF.  OF  CENTRAL  DISTANCE (USUALLY  COE=l.) 

PN-2.D0 

NP-12 

RMPG-94. / 12 . 

KSIN=30.D0*DR 
BETAP-15 . D0*DR 
COE-1 . OOODO 

C  (2)  PREDESIGN  TRANSMISSION  ERRORS  AND  CONTACT  PATH: 

C  AA-LEVEL  OF  PREDESIGNED  PARABOLIC  TRNASMISSION  ERRORS  IN  SECOND 

C  KK-COEFFICIENT  OF  CONTACT  PATH  DIRECTION  (1.D-3-CROSS  TOOTH 

C  SURFACE;  5.D6-ALONG  TOOTH  SURFACE) 

C  R-RADIUS  OF  ARC  ATTACHED  TO  CHOSEN  CONTACT  PATH 

AA-25.D0 
KK-5.D6 
R-0.3584D0 

(3)  MISALIGNMENT:  NMIS-ID  NO. (1-CROSSING  AXES,  2= INTERSECTING  AXES); 
NG-NO.  OF  MISALIGNED  ANGLES  TO  BE  SIMULATED  (FROM  -(NG~l)/2  TO 
(NG-l)/2  TIMES  GAMMAI  WITH  ODD  NG) ; 

GAMMAI- INCREMENT  OF  MISALIGNED  ANGLE  (MINUTE) ; 

NMIS-1 
NG-3 

GAMMAI -3. DO 

C  (4)  OUTPUT:  FEEI- INCREMENT  OF  ROTATION  ANGLE  OF  PINION (DEGREEE) 
FEEI-1 . 0D0*DR 
C 

C  DESCRIPTION  OF  OUTPUT  PARAMERTERS 
C 

C  FEEI -ROTATION  ANGLE  OF  PINION 

C  FEE2-ROTATION  ANGLE  OF  GEAR 

C  RP-RADIUS  OF  PINION  CONTACT  POINT 

C  RG- RADIUS  OF  GEAR  CONTACT  POINT 

C  ZP-LENGTH  OF  PINION  CONTACT  POINT  FROM  MIDDLE  SECTION 
C  ZG-LENGTH  OF  GEAR  CONTACT  POINT  FROM  MIDDLE  SECTION 
C 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 
C 

C  DEFINE  USEFUL  CONSTANTS  AND  PARAMETERS  FOR  PINION  AND  GEAR 
RMGP-l./RMPG 
NG-NP*RMPG+0 . 5 

AA-AA* (2/3600. *DR* (NP/DR/180. ) **2) 

SK-DSIN (KSIN) 

CK-DCOS (KSIN) 

SB-DSIN (BETA) 

CB-DCOS (BETA) 

TB-SB/CB 
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KSIODATAN  (SK/CK/CB) 

CKC=DCOS(KSIC) 

SKC=DSIN(KSIC) 

PT»PN*CB 

RP=NP/2./PT 

RPB-RP*CKC 

RPA-RP+1.0/PN 

RG=NG/2./PT 

RGB=RG*CKC 

RGA=RG+1.0/PN 

WRITE  (LP.56)  RP,RPB,RPA,RG,RGB,RGA,AA 
56  FORMAT  (IX, ' &&&&&& ' , 'RP=' ,F15. 7, 5X, 'RPB=' ,F15. 7, 5X,  RPA  >Fl5.7/ 
+  IX, '&&&&&& ’ , ' RG= ’ , F15 . 7, 5X, 'RGB® ' , F15 . 7 ,5X, 'RGA=  .F15.7) 

C0N1=  (1 . +SK*SK*TB*TB) 

AT  (1)  =  (CK**4/CKC**2*KK**2~SK**4*TB**2)  *RG*RG 
AT (2) = (2. *SK*SK*SB*C0Nl-2. *KK**2*SB*CK**4/CKC**2) *RG 
AT (3)  - (-2 . *CK*CK/CKC*SKC*KK**2-2 . *KK*SK*SK*TB) *RG*RG 
AT  (4)  =  (2 .  *SB*CK*CK/CKC*SKC*KK**2+2 .  *KK*CB*C0N1)  *RG 

AT(5)=CK**4/CKC*,v2*SB**2*KK,'r,f'2-CB,'c*2’'cCONl*CONl 

C*= (RP+RG) *COE 
CALL  INTMAT (S , 4 , 4) 

NGG® (NG-l)/2 
DO  505  LL=1.NG 


GAMMA=GAMMAI*FLOAT (LL-NGG) 

CG=DC0S (GAMMA/60 . *DR) 

SG=DSIN (GAMMA/60 . *DR) 

IF  (NMIS.EQ. 1)  THEN 
WRITE  (LP ,500)  COE, GAMMA 

500  FORMAT  (1H1 , ///, IX, *C-' ,F9.4, ' *(RP+RG) 

+  F8. 4, '(«)’) 

S(l,l)-CG 
S (1 , 3) *_SG 
S(3,1)®SG 
S(3,3)=CG 
ELSE 

WRITE  (LP , 501)  COE, GAMMA 

501  FORMAT  ( 1H1 , /// , IX , '  C®  1 , F4 . 2 , ' * (RP+RG) 

+  F5 . 1 ,  '  (M)  ' ) 

S(2,2)=CG 
S(2,3)— SG 
S(3,2)®SG 
S  (3, 3) ®CG 
END  IF 

NCOEF-IDINT (360 . *DR/ FEE I /FLOAT (NP) +0 . 5) 

N-2*NC0EF 

NHALF- (N+l)/2 

XIN-0. 

DO  205  L-1,2 


CROSSING  ANGLE- ' , 


INTERSECTING  ANGLE® ' , 


LSGN-(-l)**L 


DO  15  1-1, NHALF 
LI®NHALF+LSGN* (I~l) 

IF  (L.EQ.l)  NMIN-LI 
IF  (L.EQ.2)  NMAX-LI 
X  (7) -LSGN*FEEI*FLOAT (I- 1) 
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X(4)-X(7) 

X(2)»X(7)/RMPG 

X(3)»0. 

X (5)  - (- (AT (2) *X (2) +AT (4) ) +DSQRT ( (AT (2) *X (2) +AT (4) ) **2~4. *AT  (5) 
+  * (AT (3)  *X (2) +AT (1) *X (2) *X  (2) ) ) ) /2 . /AT  (5) 

X(1)«0. 

CALL  NONLIN 
X(8)=X(1)+X(2) 

C  FIND  INITIAL  VALUE  OF  X(8) 

IF  (L.EQ.l.AND.I.EQ.l)  THEN 
XIN=X(8) 

WRITE  (LP,51)  XIN 
51  FORMAT(lX, 'XIN-' .D15.7) 

GO  TO  15 
END  IF 
C 

X(8)=X(8''-XIN 
W(LI)=X  v7) /DR 
2 (LI) =X (8) /DR 

ERR (LI) - (X (8) -X (7) /RMPG) *3600 . DO/DR 
RPR(LI)=X(10) 

RGR (LI) =X  (1 1) 

2P  (Li) =X  (12) 

2G(LI)=X(13) 

C  WRITE  (LP , 20)  W(LI) ,Z(LI) ,ERR(LI) ,RPR(LI) ,RGR(LI) 

15  CONTINUE 
205  CONTINUE 

WRITE  (LP, 10) 

10  FORMAT  (////8X,  '  FEE1  (D)  '  ,8X,  ’  FEE2  (D)  '  ,  8X,  '  K- ERR  OR  (S)  ' ,  5X, 

+  1 RP ' , 13X, ' ZP ' , 1 3X , ' RG ' , 1 3X , 'ZG' ,F15.7/) 

DO  25  I*NMIN,NMAX 

WRITE  (LP,20)  W(I),Z(I),ERR(I),RPR(I),ZP(I),RGR(I),ZG(I) 

20  FORMAT  (1X.2F15. 7,F12.4,4F15. 7) 

25  CONTINUE 

NT“NMAX-NCOEF 

FII=FEEI/DR/RMPG*FLOAT(NCOEF) 

WRITE  (LP,80) 

80  FORMAT  (//,'  FIND  THE  WORKING  RANGE  FOR  ONE  TOOTH: ' , Fl5 . 7/) 

DO  55  I-NMIN.NT 

X (7)  =FEEI*FLOAT (I- (N+l) /2) /RMPG/DR 
X(8)=X(7)+FII 
KK»I+NCOEF 
ANGLE (I) =Z (KK) — Z (I) 

ERROR (I)* (ANGLE (I) -FI I) *3600. DO 
WRITE  (LP,60)  X (7) ,X (8) .ANGLE (I) .ERROR (I) 

60  FORMAT  (IX, ’ ( ' , F7 . 2 . 1 - ' f F7.2, ') : ' , F15. 7,F15. 7) 

55  CONTINUE 

DO  95  I-NMIN.NT 
ATEMP2“ERROR (I) 

IF  (I.NE.NMIN)  THEN 
IF  (ATEMP 1 *ATEMP2 . LE . 0 . DO)  GOTO  105 
END  IF 

ATEMP 1 -ATEMP 2 
95  CONTINUE 
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WRITE  (LP, 160) 

160  FORMAT  (//IX, 'MESHING  IS  DISCONTINUOUS’) 

GO  TO  505 

105  IF  (DABS (ATEMP1) .LT. DABS (ATEMP2))  1=1-1 
EMAX-0. 

EMIN-0. 

NTEMP-NCOEF+1 
DO  135  J=1 ,NTEMP 
KS-I+J-1 
ET-ERR  (KS) 

IF  (ET.LT.EMIN)  EMIN-ET 
IF  (ET.GT. EMAX)  EMAX-ET 
135  CONTINUE 

ET-EMAX-EMIN 

KK-I+NCOEF 

WRITE  (LP, 170)  Z (I) ,2 (KK) , ET 

170  FORMAT  (//IX, 'WORKING  RANGE  FOR  ONE  TOOTH:  ’,F7.2,' - *,F7.2/ 

+  IX, 'THE  MAXIMUM  KINEMATIC  ERROR:  \F12.4,’  (S) ' , 12) 

505  CONTINUE 
STOP 
END 

SUBROUTINE  FUNC 

IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,KK,KSIC 

COMMON  /BLOCKl/  X(13) ,Y(10) ,A(10, 10) ,Y1 (10) , IPVT(IO) ,WORK(10) , 
+  EPSI ,DELTA,NC,NE,NDIM 

COMMON  /BLOCK2/  S(4,4) ,AT(5) ,C,RP,RG,CK,SK,CB,SB,CKC,KSIC,TB, 

+  COEG 1 , COEG2 , RMGP , AA , DR , KK , R 

************************************************** 

20  TEMP-CKC* ( 1 . -AA*X (4) / (RMGP+1 . ) ) 

IF  (DABS (TEMP) .GT.l. DO)  THEN 
X(4)-X(7) 

GOTO  20 
END  IF 

WRITE  (6,220)  TEMP, X (4) ,AA 

220  FORMAT (IX, ' TEMP= ' ,E15 . 7 , ' X (4) * ' ,E15 . 7 , ' AA- ' , F15 . 7) 
RLAM-KSIC-DARCOS (TEMP) 

DLAMDF=-CKC*AA/ (RMGP+1 . ) /DSQRT(1 . -TEMP*TEMP) 

FEEGP=X (4) *RMGP-AA/2 . *X (4) *X (4) +RLAM 
DFGDFP=RMGP-AA*X (4) +DLAMDF 
CONA-AT (5) 

CONB-AT (2) *FEEGP+AT (4) 

CONC=AT (3) *FEEGP+AT(1) *FEEGP*FEEGP 
COND-CONB*CONB-4 . *CONA*CONC 
IF  (COND.LT.O.)  THEN 
X(4)«X(7) 

GOTO  20 
END  IF 

COND-DSQRT (COND) 

TP- (-CONB+COND) /2 . /CONA 
IF  (COND. EQ. 0.0)  THEN 
DTPDFG— AT  (2)  /2 .  /CONA 
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ELSE 

DTPDFG-  (-AT (2)  -  (AT  (2)  *C0NB-2. *CONA*  (AT (3)  +2 . *AT (1)  *FEEGP) )  /COND) 
+  /2./CONA 

END  IF 

£  *******  £*:&************&  *********** 

C  X(4) =FEEP;  X  (7) =FEE1 
AP1=X(4)-X(7) 

AP2=AP1-RLAM 
AP3=AP2+KSIC 
AP4=X (4) -RLAM+KSIC 
DAP4DF= 1 . -DLAMDF 
CAP4=DCOS (AP4) 

SAP4=DSIN(AP4) 

AP—TP*SB+RG*FEEGP 

DAPDF*=  (-SB*DTPDFG+RG)  *DFGDFP 

MU=X (4) +KSIC-RLAM 

DMUDFP= 1 . +DLAMDF 

ALP1=R* (DSIN (MU+X (3) ) -DSIN (MU) ) 

ALP2=R* (DCOS (MU+X (3) ) -DCOS (MU) ) 

XPF=CK*CK/CKC*DCOS (AP3) *AP+RG*DSIN (AP2) - (RP+RG) *DSIN (API) +ALP2 
YPF=CK*CK/CKC*DSIN (AP3) *AP-RG*DCOS (AP2) + (RP+RG) *DCOS (API) +ALP1 
ZPF=TP*CB*CK*CK/CKC/CKC-SK*SK*TB*RG*FEEGP 
c********************************************************** 

DXDA=-DSIN (MU+X  (3) ) 

DYDA=DCOS (MU+X  (3) ) 

DXDF=CK*CK/CKC* (CAP4*DAPDF"SAP4*AP*DAP4DF) +RG*DCOS (X (4) -RLAM) 

+  * (1 . -DLAMDF) - (RP+RG) *DCOS (X (4) ) -ALP1*DMUDFP 

DYDF=CK*CK/CKC* (SAP4*DAPDF+CAP4*AP*DAP4DF) +RG*DSIN (X (4) -RLAM) 

+  *(1. -DLAMDF) - (RP+RG) *DSIN (X (4) ) +ALP2*DMUDFP 

DZDF=CB*CK*CK/CKC/CKC*DTPDFG*DFGDFP-SK*SK*TB*RG*DFGDFP 
XNP=DYDA*DZDF 
YNP— DXDA*DZDF 
ZNP=DXDA*DYDF-DYDA*DXDF 
RMN=DSQRT(XNP*XNP+YNP*YNP+ZNP*ZNP) 

C  WRITE (6 , 130)  DXDA , DZDA , DXDF , DYDF , DZDF , XNP , YNP , ZNP , RMN , TP , AP 
C  130  FORMAT (IX, ' $$$ ' ,5F15.7/4X,5F15.7) 

XNP=XNP/RMN 

YNP-YNP/RMN 

ZNP=ZNP/RMN 

XNPF=XNP*DCOS (X (7) ) +YNP*DSIN (X  (7) ) 

YNPF*YNP*DCOS (X (7) ) ~XNP*DS IN (X (7) ) 

ZNPF=ZNP 

CF2FG“DCOS(X(l)) 

SF2FG»DSIN (X (1) ) 

CF2FGK-DCOS (X (1) +KSIC) 

SF2FGK-DS IN (X ( 1) +KS IC) 

AG— X(5)*SB+RG*X(2) 

RFG1-CK*CK/CKC*CF2FGK*AG+RG*SF2FG 

RFG2-CK*CK/CKC*SF2FGK*AG-RG*CF2FG 

RFG3-X  (5) *CB-AG*SK*SK*TB 

RNFG1-CK*CB*CF2FG-SK*SF2FG 

RNFG2-CK*CB*SF2FG+SK*CF2FG 

RNFG3"CK*SB 
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T 


XGF=RFG1*S (1 , 1) +RFG2*S (1,2) +RFG3*S  (1,3) 

YGF-RFG1*S (2 , 1) +RFG2*S (2 , 2) +RFG3*S  (2,3) 

ZGF=RFG1*S (3,1) +RFG2*S (3 , 2) +RFG3*S (3 , 3) 

XNGF-RNFG1*S (1 , 1) +RNFG2*S (1,2) +RNFG3*S (1 , 3) 

YNGF=RNFG1*S (2,1) +RNFG2*S (2 , 2) +RNFG3*S (2 , 3) 

ZNGF-RNFG1*S (3,1) +RNFG2*S (3 , 2) +RNFG3*S  (3 , 3) 

WRITE  (6,100)  X(l)  ,X(2)  ,X(3)  ,X(4)  ,X(5) 

WRITE  (6,100)  XPF , YPF , ZPF , XGF , YGF , ZGF 
100  FORMAT  (IX,  'imr  ,8E15.7) 

WRITE  (6,100)  XNPF , YNPF , ZNPF , XNGF , YNGF , ZNGF 

Y (1) “XPF-XGF 

Y(2)=YPF-YGF-C 

Y(3)=ZPF-ZGF 

Y (5) “XNPF-XNGF 

Y(4)=ZNPF-ZNGF 

X ( 10) =DSQRT (XPF*XPF+YPF*YPF) 

X(ll) =DSQRT (XGF*XGF+YGF*YGF) 

X(12)=ZPF 

X(13)*=ZGF 

WRITE  (6,20)  (Y  (II) ,11  =  1,5) 

20  FORMAT  (IX, '$$$$' .6F15.7) 

RETURN 

END 

SUBROUTINE  INTMAT  (A , N , M) 

THIS  SUBROUTINE  IS  USED  TO  INITIATE  THE  MATRIX,  WITH  UNIT  DIAGONAL 
ELEMENTS  AND  NULL  OTHER  ELEMENTS 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(4,4) 

DO  5  1=1, N 
DO  5  J=1,M 
A(I,J)-0. 

IF  (I.EQ.J)  A(I , j) =1 . 

5  CONTINUE 
RETURN 
END 


SUBROUTINE  NONLIN 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /BLOCK1/  X(13) ,Y(10) ,A(10, 10) ,Y1 (10) , IPVT(IO) ,WORK(10) , 
+  EPSI, DELTA, NC.NE.NDIM 

C 

DO  5  1=1, NC 
CALL  FUNC 

C  WRITE  (6,10)  I,  (X(J) ,Y(J) , J-1,5) 

C  10  FORMAT (IX, ***** , 15/5 (IX, 2D15 . 7/) ) 

DO  15  J=!,NE 

IF  (DABS (Y (J) ) .GT. EPSI)  GO  TO  25 
15  CONTINUE 
GO  TO  105 
25  DO  35  J=1,NE 
35  Y1(J)=Y(J) 


r  -V 
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DO  45  J=1,NE 
DIFF-DABS (X (J) ) *DELTA 
IF  (DABS(X(J)) .LT. l.D-12)  DIFF-DELTA 
XMAM-X(J) 

X(J)«X(J)-DIFF 
CALL  FUNC 
X  (  J)  *=XMAM 
DO  55  K-l.NE 

55  A(K, J) “ (Y1 (K) -Y (K) ) /DIFF 
45  CONTINUE 
DO  65  J=1 ,NE 
65  Y(J)—  Y1  (J) 

C  DO  205  K=1 ,NE 
C  205  WRITE  (6,245)  (A(K, J) , J=1 ,NE) 

C  245  FORMAT  (IX,  5D15 . 7) 

CALL  DECOMP  (NDIM.NE.A.COND, IP VT, WORK) 
CALL  SOLVE  (NDIM,NE,A, Y, IPVT) 

DO  75  J=1 ,NE 
75  X  (J)  ”X  (J)  +Y  (J) 

5  CONTINUE 
C  105  WRITE  (6,20)  I 
C  20  FORMAT  (IX,' I- ’,12) 

105  RETURN 
END 


SUBROUTINE  DECOMP  (NDIM.N, A.COND, IPVT, WORK) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (NDIM.N) , WORK (N) , IPVT (N) 

C 

C  DECOMPOSES  AREAL  MATRIX  BY  GAUSSIAN  ELIMINATION, 

C  AND  ESTIMATES  THE  CONDITION  OF  THE  MATRIX. 

C 

C  -COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
C  M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL,  1977) 

C 

C  USE  SUBROUTINE  SOLVE  TO  COMPUTE  SOLUTIONS  TO  LINEAR  SYSTEM. 

C 

C  INPUT. . 

C 

C  NDIM  -  DECLARED  ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  A 

C  N  ORDER  OF  THE  MATRIX 

C  A  MATRIX  TO  BE  TRIANGULAR I ZED 

C 

C  OUTPUT.. 

C 

C  A  CONTAINS  AN  UPPER  TRIANGULAR  MATRIX  U  AND  A  PREMUTED 
C  VERSION  OF  A  LOWER  TRIANGULAR  MATRIX  I-L  SO  THAT 

C  (PERMUTATION  MATRIX) * A- L*U 

C 

C  COND  -  AN  ESTIMATE  OF  THE  CONDITION  OF  A. 

C  FOR  THE  LINEAR  SYSTEM  A*X  -  B  ,  CHANGES  IN  A  AND  B 
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C  MAY  CAUSE  CHANGES  COND  TIMES  AS  LARGE  IN  X. 

C  IF  COND+l.O  .EQ.  COND  ,  A  IS  SINGULAR  TO  WORKING 

C  PRECISION.  COND  IS  SET  TO  1.0D+32  IF  EXACT 

C  SINGULARITY  IS  DETECTED. 

C 

C  IPVT  -  THE  PIVOT  VECTOR 

C  IPVT (K)  -  THE  INDEX  OF  THE  K-TH  PIVOT  ROW 

C  IPVT (N)  -  (-1)** (NUMBER  OF  INTERCHANGES) 

C 

C  WORK  SPACE..  THE  VECTOR  WORK  MUST  BE  DECLARED  AND  INCLUDED 
C  IN  THE  CALL.  ITS  INPUT  CONTENTS  ARE  IGNORED. 

C  ITS  OUTPUT  CONTENTS  ARE  USUALLY  UNIMPORTANT. 

C 

C  THE  DETERMINANT  OF  A  CAN  BE  OBTAINED  ON  OUTPUT  BY 
C  DET(A)  -  IPVT (N)  *  A(1 , 1)  *  A(2,2)  *  ...  *  A(N,N)  . 

C 

IPVT (N) “1 

IF  (N.EQ.l)  GO  TO  150 
NM1=N-1 

C  COMPUTE  THE  1-NORM  OF  A  . 

ANORM=0 .  DO 
DO  20  J=1 ,N 
T=0.D0 
DO  10  1=1, N 
10  T“T+DABS (A(I , J) ) 

IF  (T.GT.ANORM)  ANORM=T 
20  CONTINUE 

DO  GAUSSIAN  ELIMINATION  WITH  PARTIAL 
PIVOTING. 

DO  70  K=1,NM1 
KP1-K+1 

FIND  THE  PIVOT. 

M=K 

DO  30  I=KP1,N 

IF  (DABS (A(I ,K) ) . GT. DABS (A(M,K) ) )  M=I 
30  CONTINUE 
IPVT  (K)  *=M 

IF  (M.NE.K)  IPVT (N) “-IPVT (N) 

t=a(m,k) 
a(m,k) -a(k,k) 
a(k,k)=t 

SKIP  THE  ELIMINATION  STEP  IF  PIVOT  IS  ZERO. 
IF  (T.EQ.O.DO)  GO  TO  70 

COMPUTE  THE  MULTIPLIERS. 


INTERCHANGE  AND  ELIMINATE  BY  COLUMNS. 


GO  TO  60 
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DO  40  I-KPl.N 
40  A(I,K)—  A(I,K)/T 

C 

DO  60  J-KPl.N 
T-A(M, J) 
A(M,J)-A(K,J) 
A(K, J)-T 
IF  (T.EQ.O.DO) 
DO  50  I-KPl.N 


n  o  non 


50  A(I,J)-A(I,J)+A(I,K)*T 

60  CONTINUE 
70  CONTINUE 
C 

C  COND  -  (1-NORM  OF  A)* (AN  ESTIMATE  OF  THE  1-NORM  OF  A-INVERSE) 

C  THE  ESTIMATE  IS  OBTAINED  BY  ONE  STEP  OF  INVERSE  ITERATION  FOR  THE 
C  SMALL  SINGULAR  VECTOR.  THIS  INVOLVES  SOLVING  TWO  SYSTEMS 
C  OF  EQUATIONS,  (A-TRANSPOSE) *Y  -  E  AND  A*Z  =  Y  WHERE  E 
C  IS  A  VECTOR  OF  +1  OR  -1  COMPONENTS  CHOSEN  TO  CAUSS  GROWTH  IN  Y. 

C  ESTIMATE  -  (1-NORM  OF  Z) / (1-NORM  OF  Y) 

C 

C  SOLVE  (A-TRANSPOSE) *Y  “  E  . 

DO  100  K=1,N 
T=0.D0 

IF  (K.EQ.l)  GO  TO  90 
KM1-K-1 
DO  80  1=1, KM1 
80  T=T+A(l,K)*WORK(l) 

90  EK= 1 . DO 

IF  (T.LT.O.DO)  EK=-1 .DO 
IF  (A(K,K) .EQ.O.DO)  GO  TO  160 
100  WORK (K) =- (EK+T) /A (K , K) 

DO  120  KB-l.NMl 
K=N-KB 
T=0.D0 
KP1-K+1 

DO  110  I=KP 1 , N 
110  T=T+A(l,K)*WORK(K) 

WORK (K) =T 
M=IPVT(K) 

IF  (M.EQ.K)  GO  TO  120 
T=WORK(M) 

WORK (M) -WORK (K) 

WORK (K) “T 
120  CONTINUE 
C 

YNORM-O.DO 
DO  130  1=1, N 

130  YNORM- YNORM+DAB S(WORK(l)) 

SOLVE  A*Z  =  Y 

CALL  SOLVE  (NDIM , N , A , WORK , I P VT) 

ZNORM-0 . DO 
DO  140  1=1, N 

140  ZNORM=ZNORM+DABS (WORK(l) ) 

ESTIMATE  THE  CONDITION. 

COND-ANORM*ZNORM/ YNORM 
IF  (COND.LT. 1.D0)  COND- 1. DO 
RETURN 

C  1-BY-l  CASE.. 

150  COND- 1. DO 

IF  (A(l , 1) . NE.O.DO)  RETURN 
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EXACT  SINGULARITY 


160  COND-1.0D32 
RETURN 
END 

SUBROUTINE  SOLVE  (NDIM,N, A, B , IPVT) 


IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  A(NDIM,N) ,B(N) ,IPVT(N) 


SOLVES  A  LINEAR  SYSTEM,  A*X  =  B 

DO  NOT  SOLVE  THE  SYSTEM  IF  DECOMP  HAS  DETECTED  SINGULARITY. 


-COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL,  1977) 


INPUT. . 

NDIM  =  DECLARED  ROW  DIMENSION  OF  ARRAY  CONTAINING  A 
N  -  ORDER  OF  MATRIX 

A  =  TR I ANGULAR IZED  MATRIX  OBTAINED  FROM  SUBROUTINE  DECOMP 

B  =  RIGHT  HAND  SIDE  VECTOR 

IPVT  =  PIVOT  VECTOR  OBTAINED  FROM  DECOMP 

OUTPUT. . 

B  =  SOLUTION  VECTOR,  X 


DO  THE  FORWARD  ELIMINATION. 

IF  (N.EQ.l)  GO  TO  50 
NM1°N-1 
DO  20  K= 1 , NM 1 
KP 1 =K+ 1 
M=IPVT(K) 

T=B (M) 

B(M)«B  (K) 

B (K) «T 

DO  10  I=KP1,N 
10  B(I)-B(I)+A(I,K)*T 

20  CONTINUE 

NOW  DO  THE  BACK  SUBSTITUTION. 

DO  40  KB=1,NM1 
KM1-N-KB 
K-KM1+1 

B(K)=B  (K)  /A(K,K) 

T— B  (K) 

DO  30  I-l.KMl 
30  B(I)-B(I)+A(I,K)*T 

40  CONTINUE 
50  B(l)-B(l)/A(l,l) 

RETURN 

END 


135 


FLOWCHART  FOR  PROGRAM  VI 


l 


t 


t 


f 

i 


136 


n  n  o  n  n 


**************************** *************************** 


c. . 
c. . 
c. . 
c. . 
c. . 
c. . 
c. . 
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c. . 
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c. . 
c 


*  * 

*  * 

*  PROGRAM  VI  * 

*  TRANSMISSION  ERRORS  OF  HELICAL  GEARS  WITH  THEIR  * 

*  AXES  DEFORMED  BY  INTERACTING  FORCE  * 

*  * 

*  AUTHORS:  FAYDOR  LITVIN  * 

*  JIAO  ZHANG  * 

*  * 

*  * 


******************************************************* 


C  PURPOSE 


C 


C  THIS  PROGRAM  IS  USED  TO  CALCULATE  THE  TRANSMISSION  ERRORS  OF  A 
C  HELICAL  PINION  AND  A  HELICAL  GEAR  IN  MESHING  WHEN  THEIR  AXES  ARE 
C  DEFORMED  BY  INTERACTING  FORCE (BOTH  PINION  AND  GEAR  ARE  NOT  CROWNED) 
C 

C  NOTE 


C 

C  THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  IT  CAN  BE  COMPILED  BY  V 
C  COMPILER  IN  IBM  MAINFRAME  OR  FORTRAN  COMPILER  IN  VAX  SYSTEM. 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,MUO,KSlC 

DIMENSION  Z (99) .ANGLE (99) .ERROR (99) .ERR (99) ,W(99) ,RPR(99) ,RGR(99)  , 
+  S 1  (A ,  4)  ,  S2  (4 , 4)  ,  ZP  (99) 

COMMON  /B LOCK1/  X  (12) , Y (10) , A (10, 10) , Y1 (10) , IPVT (10) .WORK(IO) , 

+  EPSI, DELTA, NC,NE,NDIM,NCTL,CX2 

COMMON  /BLOCK2/  S(4,4) .C.RP.RG.CK.SK.CB.SB.CKC.KSIC.ZG.TB.ZNPF, 

+  COEG1 .COEG2.CB1 , SB1 , TB1 


DEFINE  PARAMETERS  USED  BY  PROGRAMS 


(I)  IN  ANG  LP  ARE  UNIT  NUMBERS  ASSIGNED  TO  THE  INPUT  AND  OUTPUT 
DEVICES 

IN-5 
LP-6 

C  (2)  NDBUG  IS  USE  TO  CONTROL  THE  AUXILIARY  OUTPUT  FOR  DEBUGGING 
NDBUG-1 

C  (3)  NC  IS  THE  UPPER  LIMITATION  OF  REPEATATION  FOR  SOLVING  NONLINEAR 
C  EQUTIONS; 

C  EPSI  IS  THE  CLEARANCE  OF  FUNCTION  VALUES  WHEN  THE  FUNCTIONS 

C  IS  CONSIDERED  AS  SOLVED  (ALL  FUNTIONS  HAVE  FORMS  OF  F(X)-O); 

C  DELTA  IS  THE  RELATIVE  DIFFERENCE  FOR  TAKING  DERIVATIVES 

C  NC,  EPSI  AND  DELTA  MAY  BE  CHANGED  WHEN  SOLUTIONS  ARE  DIVERGENT 

C  OR  LESS  ACCURATE 

NC-100 
DELTA- l.D-4 
EPSI-l.D-13 

C  (4)  OTHER  PARAMETERS (DON ' T  CHANGE) 

NDIM-10 

DR-DATAN ( 1 . DO) /45 . DO 
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C  DEFINE  INPUT  PARAMTERS  OF  PROBLEM (USE  INCH  AS  UNIT  OF  LENGTH) 

C  (1)  PINION  AND  GEAR:  PN=DIAMETRAL  PITCH;  NP=PINION  TOOTH  NUMBER; 
C  RMPG=TOOTH  NUMBER  RAT 10 (GEAR  TOOTH  NO./NP); 

C  KSIN=PRESSURE  ANGLE  IN  NORMAL  SECTION; 

C  BETAP=HELIX  ANGLE  OF  PINION  AND  GEAR; 

C  PAD= HEIGHT  OF  ADDENDUM  OF  PINION; 

C  GAD=HEIGHT  OF  ADDENDUM  OF  GEAR; 

C  ZG=GEAR  TOOTH  LENGTH  (PINION  TOOTH  LENGTH  IS  LONGER) 

C  COE=COEFF.  OF  CENTRAL  DISTANCE (USUALLY  C0E=1.) 

PN=10.D0 
NP=20 
RMPG=2.DC 
KSIN=20.D0*DR 
BETAP=15.D0*DR 
PAD=1 . O/PN 
GAD=1 . O/PN 
2G=  5./PN 
C0E=1 . OOODO 

C  (2)  SHAFT  DEFORMATION: 

C  NSIM=MOD£L  ID  NO. (1=SIMPLIFIED  DEFORMATION  MATRIX; 

C  2=UNSIMPLIFIED  DEFORMATION  MATRIX) 

C  RLAMP=PINION  SHAFT  ROTATION 

C  RLAMG=GEAR  SHAFT  ROTATION 

C  RVP=PINION  SHAFT  DEFLECTION 

C  RVG=GEAR  SHAFT  DEFLECTION 

C 

NSIM=1 

RLAMP=  2./60.*DR 
RLAMG-  2./60.*DR 
RVP=0.0125 
RVG=0.0125 

C  (3)  OUTPUT:  FEEI= INCREMENT  OF  ROTATION  ANGLE  OF  PINION (DEGREEE) 
FEE 1=1 . ODO*DR 
C 

C  DESCRIPTION  OF  OUTPUT  PARAMERTERS 
C 

C  FEE1-R0TATI0N  ANGLE  OF  PINION 

C  FEE2=R0TATI0N  ANGLE  OF  GEAR 

C  RP=RADIUS  OF  PINION  CONTACT  POINT 

C  RG= RADIUS  OF  GEAR  CONTACT  POINT 
C 

C  FIND  AUXILIARY  VALUES  FOR  CALCULATION 
C 

c 

C  DEFINE  USEFUL  CONSTANTS  AND  PARAMETERS  FOR  PINION  AND  GEAR 
DELTAB«DR/60.*(  O.DO) 

BETAP 1 -BETAP+DELTAB 
CLP-DCOS(RLAMP) 

SLP=DSIN (RLAMP) 

CLG-DCOS (RLAMG) 

SLG-DS IN (RLAMG) 

WRITE  (LP,4)  CLP,SLP,CLG,SLG 
4  FORMAT (IX, AF15 . 7) 

SK-DSIN(KSIN) 
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CK-DCOS(KSIN) 

SB-DS IN (BETAP) 

CB=DCOS (BETAP) 

TB-SB/CB 
SB1=DSIN (BETAP1) 

CBl-DCOS (BETAP 1) 

TB1-SB1/CB1 
KSIC=DATAN (SK/CK/CB) 

CKC=DCOS (KSIC) 

SKC=DSIN (KSIC) 

PT«PN*CB 
RP=NP/2./PT 
RPA=RP+PAD 
RPATOL=RPA-0. 0005D0 
RG=RP*RMPG 
RGA=RG+GAD 
RGATOL=RGA-0 . 0005D0 
WRITE  (LP,56)  RPATOL.RGATOL 
56  FORMAT  (IX, ' &&&&&& 1 , ' RPA= ' , F15 . 7 , 5X, ' RGA= ' , F15 . 7) 

COEGl=l . +SK*SK*TB*TB 
C0EG2=RG/C0EG1 
ZNPF=CK*SB1 
C= (RP+RG) *COE 
CALL  INTMAT (S , 4, 4) 

CALL  INTMAT (SI, A, 4) 

CALL  INTMAT (S2, A, A) 

DO  505  LL=1 ,2 
NSIM=LL 

IF  (NSIM.EQ. 1)  THEN 
WRITE  (LP.500) 

500  FORMAT  (1H1,///,1X, 

+  'THE  CASE  OF  SIMPLIFIED  DEFORMATION  MATRIX  OF  GEAR  AXES') 

S ( 1 , 3) - (RLAMP+RLAMG) *CKC 
S (2 , 3) - (RLAMP+RLAMG) *SKC 
S  (3, 1)  =-S  (1 , 3) 

S (3 , 2)  =-S (2 , 3) 

S ( 1 ,  A)  =  (R VP+R VG) *CKC 
S (2 , A) - (RVP+RVG) *SKC+C 
S (3 , A) — C*RLAMP*SKC 
ELSE 

WRITE  (LP.501) 

501  FORMAT  (1H1,///,1X, 

+  'THE  CASE  OF  UNSIMPLIFIED  DEFORMATION  MATRIX  OF  GEAR  AXES') 

SI (1 , 1) =1 . +CKC*CKC* (CLP-1 .) 

S 1 ( 1 , 2) -SKC*CKC*  (CLP-1 . ) 

S1(1,3)-CKC*SLP 
SI (2, 1)*S1 (1,2) 

S 1 (2 , 2) - 1 . +SKC*SKC* (CLP- 1 . ) 

S 1  (2 , 3) “SKC*SLP 
SI  (3,1)—  Sl(l, 3) 

SI  (3, 2) —SI  (2,3) 

SI (3,3) “CLP 

Sl(l, A) -RVP*CKC*CLP 

S1(2,4)-RVP*SKC*CLP 
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51  (3,4)— RVP*SLP 

S2 (1 , 1) - 1 . +CKC*CKC* (CLG- 1 . ) 

52  (1 , 2)  -SKC*CKC*  (C.LG-1 . ) 

S2 (1 , 3) =CKC*SLG 

S2 (2, 1) “SI  (1 , 2) 

S2(2,2)-l. +SKC*SKC* (CLG-1 . ) 

S2 (2 , 3) =SKC*SLG 

S2(3,l)—  Sl(l, 3) 

S2(3 , 2) —SI (2,3) 

S2  (3 , 3) “CLG 
S2  (1 , 4) =RVG*CKC 
S2  (2 , 4) =RVG*SKC+C 
S2(3,4)-0. 

DO  12  MMI  =  1 ,3 
DO  12  MMJ-1,4 
S (MMI ,MMJ) “0. 

DO  12  MMK=1,4 

S  (MMI ,MMJ) =S (MMI ,MMJ) +S1 (MMI ,MMK) *S2 (MMK.MMJ) 

12  CONTINUE 
END  IF 

NCOEF* IDINT (360 . *DR/FEEI /FLOAT (NP) +0 . 5) 

N=3*NCOEF 
NHALF- (N+l)/2 
XIN=0. 

DO  205  L-1,2 
LSGN= (-1)**L 
NCTL=0 
NE=4 

DO  5  1=1 ,NE 
5  X(I)=0.D0 
X(10)=0. 

X(11)=0. 

DO  15  1=1 , NHALF 
IF  (NCTL.EQ.2)  GOTO  205 
LI-NHALF+LSGN* ( I — 1 ) 

IF  (L.EQ.l)  NMIN=LI 
IF  (L.EQ.2)  NMAX=LI 
X(7)=LSGN*FEEI*FLOAT(I-l) 

IF  (X(10) .LE.RPATOL)  THEN 
X(4)=- X(7) 

ELSE  "tv 

NE=3 

NCTL-NCTL+1 
END  IF 

IF  (X(ll) .LT.RGATOL)  THEN 
X(2)«X(7) /RMPG 
ELSE 

NCTL-NCTL+1 

CX2-X(2) 

END  IF 

X (3) - (ZG-RP*X (4) *SK*SK*TB) /CB /COEG1 

CALL  NONLIN 

X(8)-X(l)+X(2) 

C  FIND  INITIAL  VALUE  OF  X(8) 
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IF  (L.EQ.l.AND.I.EQ.l)  THEN 
XIN=X  (8) 

GO  TO  15 
END  IF 

X(8)*=X(8)-XIN 
W(LI)-X(7)/DR 
Z (LI) -X (8) /DR 

ERR  (LI)  =*  (X  (8)  -X  (7)  /RMPG)  *3600 .  DO/DR 
RPR  (LI)  =X  (10) 

RGR(LI)-X(11) 

ZP (LI) =X (12) 

WRITE  (LP,21)  LI,W(LI) ,Z(LI) , ERR (LI) ,RPR(LI) ,RGR(LI) 

21  FORMAT  (1X,I2,1X,5F15.7,F15.7) 

15  CONTINUE 
205  CONTINUE 

WRITE  (LP , 10) 

10  FORMAT  (////8X, 'FEEl(D) ' ,8X, ' FEE2 (D) ’ ,8X, 'K-ERROR(S) ' ,5X, 

+  'RP' , 13X, 'RG' ,F15.7/) 

DO  25  I=NMIN,NMAX 

WRITE  (LP , 20)  W(I)  ,Z(I)  ,ERR(I)  ,RPR(I)  ,RGR(I) 

20  FORMAT  (IX, 5F15 . 7 , F15 . 7) 

25  CONTINUE 

NT»NMAX-NCOEF 

FII=FEEI/DR/2 . DO*FLOAT (NCOEF) 

WRITE  (LP ,80) 

80  FORMAT  (//,'  FIND  THE  WORKING  RANGE  FOR  ONE  TOOTH: ' , F15 . 7/) 
DO  55  I=NMIN,NT 

X  (7) -FEEI*FLOAT ( I- (N+ 1) /2) /DR/2 . DO 

X(8)=X(7)+FII 

KK-I+NCOEF 

ANGLE (I) “Z (KK) ~Z (I) 

ERROR (I) - (ANGLE (I) -FI I) *3600 . DO 

WRITE  (LP ,60)  X(7) ,X(8) ,ANGLE(I) ,ERROR(l) 

60  FORMAT  (IX, '  (' ,F7.2, ' - ' , F7. 2, ' ) : ' , F15. 7 , F15 . 7) 

55  CONTINUE 

DO  95  I=NMIN,NT 
ATEMP2-ERROR (I) 

IF  (I.NE.NMIN)  THEN 

IF  (ATEMP 1 *ATEMP2 . LE • 0 . DO)  GOTO  105 

END  IF 

ATEMP 1“ ATEMP 2 
95  CONTINUE 

WRITE  (LP, 160) 

160  FORMAT  (//IX, 'MESHING  IS  DISCONTINUOUS') 

GO  TO  505 

105  IF  (DABS (ATEMP1) . LT. DABS (ATEMP2) )  I-I-l 
EMAX-0. 

EMIN-0. 

NTEMP-NCOEF+1 
DO  135  J-l.NTEMP 
KS-I+J-1 
ET-ERR (KS) 

IF  (ET.LT.EMIN)  EMIN-ET 


141 


IF  (ET.GT.EMAX)  EMAX-ET 
135  CONTINUE 

ET-EMAX-EMIN 

KK-I+NCOEF 

WRITE  (LP, 170)  Z(I) ,Z(KK) ,ET 

170  FORMAT  (//IX, 'WORKING  RANGE  FOR  ONE  TOOTH:  ' ,F7.2,' - ' ,F7.2/ 

+  IX, 'THE  MAXIMUM  KINEMATIC  ERROR:  ',F15.7,'  (S) ' , 12) 

505  CONTINUE 
STOP 
END 
C 

SUBROUTINE  FUNC 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLE  PRECISION  KSIN,MU,MUO,KSIC 

COMMON  /BLOCK1/  X (12) , Y (10) ,A(10, 10) , Y1 (10) , IPVT (10) , WORK(IO) , 

+  EPSI, DELTA, NC,NE,NDIM,NCTL,CX2 

COMMON  /BLOCK2/  S (4, 4) , C , RP ,RG , CK, SK, CB , SB , CKC, KSIC.ZG, TB ,ZNPF, 
+  COEG1 ,COEG2 , CB1 , SB1 , TBl 

C  X(4)=FEEP;  X (7) =FEE1 
X(6)=X(4)+X(7) 

CF1FP=DC0S (X(6) ) 

SF1FP=DSIN (X  (6) ) 

AP=X (3) *SB 1+RP*X (4) 

XPF=-CK*CK/CKC*DCOS  (X(6)  -KSIC)  *AP+RP*SF1FP 
YPF=CK*CK/CKC*DSIN (X (6) -KSIC) *AP+RP*CF1FP 
ZPF=X(3)*CB1+AP*SK*SK*TB1 
XNPF=CK*CB 1 *CF 1 FP+SK*SF 1 FP 
YNPF=-CK*CB1*SF1FP+SK*CF1FP 
C  ZNPF=CK*SB1 

CF2FG=DC0S (X (1)  ) 

SF2FG=DSIN(X(1)) 

CF2FGK=DC0S (X (1)  +KSIC) 

SF2FGK=DSIN (X (1) +KSIC) 

AG1=  (-ZG*TB+RG’VX  (2)  )  /COEGl 
RFG1=CK,VCK/CKC*CF2FGK,'AG1+RG,'CSF2FG 
RFG2=CK*CK/CKC*SF2FGK*AG1-RG*CF2FG 
C  RFG3=ZG 

RTFG1=CK*CK/CKC* (  AG1*SF2FGK+C0EG2*CF2FGK) -RG*CF2FG 
RTFG2=CK*CK/CKC* (-AG1*CF2FGK+C0EG2*SF2FGK) -RG*SF2FG 
C  RTFG3=0. 

XGF=RFG1*S (1 , 1) +RFG2*S (1 , 2) +ZG*S (1 , 3) +S (1 , 4) 

YGF=RFG1*S (2 , 1) +RFG2*S (2 , 2) +ZG*S (2 , 3) +S (2 , 4) 

ZGF=RFJ1*S  (3 , 1)  +RFG2y'S  (3 , 2)  +ZG*S  (3 , 3)  +S  (3 , 4) 

XTGF=RTFG1*S (1,1) +RTFG2*S (1,2) 

YTGF=RTFG1*S (2,1) +RTFG2*S  (2,2) 

ZTGF=RTFG1*S (3,1) +RTFG2*S  (3,2) 

C  WRITE  (6,100)  RTFG1 , RTFG2 , AG1 ,COEG2, COEGl , CF2FGK, SF2FGK 
C  WRITE  (6,100)  X (1) ,X(2) ,X(3) ,X(4) ,X(5) 

C  WRITE  (6,100)  XPF , YPF , ZPF , XGF , YGF , ZGF 

C  100  FORMAT  (IX,  "LVLVL'  .8E15.7) 

C  WRITE  (6,100)  XNPF , YNPF , ZNPF , XTGF , YTGF , ZTGF 
Y(l) “XPF-XGF 
Y  (2) =YPF-YGF 
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Y (3) “ZPF-ZGF 
IF  (NCTL.NE.O)  THEN 
Y  (4)  «X(2)-CX2 
ELSE 

Y (4) =XNPF*XTGF+ YNPF*YTGF+ZNPF*ZTGF 
END  IF 

X (10) =DSQRT (XPF*XPF+YPF*YPF) 

X  (1 1)  =DSQRT (RFG1*RFG1+RFG2*RFG2) 

X  (12) =ZPF 

WRITE  (6,20)  (Y(II) , 11=1,4) 

20  FORMAT  (IX, ’$$$$' ,6F15.7) 

RETURN 
END 

SUBROUTINE  INTMAT  (A,N,M) 

THIS  SUBROUTINE  IS  USED  TO  INITIATE  THE  MATRIX,  WITH  UNIT  DIAGONAL 
ELEMENTS  AND  NULL  OTHER  ELEMENTS 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (4, 4) 

DO  5  1=1,  N 
DO  5  J=1,M 
A(I,  J)=0. 

IF  (I.EQ.J)  A(I , J) =1 . 

5  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  NONLIN 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /BLOCK1/  X  (12) , Y  (10) , A (10, 10)  , Yl (10) , IPVT (10) , WORK ( 10) , 

+  EPS I , DELTA, NC,NE,NDIM,NCTL,CX2 

C 

DO  5  1=1, NC 
CALL  FUNC 

C  WRITE  (6,10)  I, (X(J) ,Y(J) , J=l,4) 

C  10  FORMAT (IX, '***', 15/5 (IX, 2D15. 7/)) 

DO  15  J=1 ,NE 

IF  (DABS (Y (J) ) .GT.EPSI)  GO  TO  25 
15  CONTINUE 
GO  TO  105 
25  DO  35  J=1 ,NE 
35  Yl (J) =Y (j) 

DO  45  J=1 ,NE 
DIFF=DABS (X ( J) ) *DELTA 
IF  (DABS(X(J)) .LT. l.D-12)  DIFF=DELTA 
XMAM=X ( J) 

X(J)=X(J)-DIFF 
CALL  FUNC 
X (J) =XMAM 
DO  55  K=1 ,NE 

55  A(K,J)=(Y1(K)-Y(K))/DIFF 
45  CONTINUE 
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DO  65  J-l.NE 
65  Y  (J) “-Y1 (J) 

C  DO  85  K-l.NE 

C  85  WRITE  (6,104)  (A (K, J) , J-l.NE) ,Y(K) 

C  104  FORMAT (IX, 'A' , 2X, 5E15 . 7) 

CALL  DECOMP  (NDIM.NE.A.COND, IPVT.WORK) 
CALL  SOLVE  (NDIM.NE, A, Y , IPVT) 

DO  75  J-l.NE 
75  X(J)-X(J)+Y(J) 

5  CONTINUE 
105  RETURN 
END 


SUBROUTINE  DECOMP  (NDIM, N, A, COND, IPVT.WORK) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (NDIM, N) , WORK (N) , IPVT (N) 

C 

C  DECOMPOSES  AREAL  MATRIX  BY  GAUSSIAN  ELIMINATION, 

C  AND  ESTIMATES  THE  CONDITION  OF  THE  MATRIX. 

C 

C  -COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
C  M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL,  1977) 

C 

C  USE  SUBROUTINE  SOLVE  TO  COMPUTE  SOLUTIONS  TO  LINEAR  SYSTEM. 

C 

C  INPUT.. 

C 

C  NDIM  -  DECLARED  ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  A 

C  N  ORDER  OF  THE  MATRIX 

C  A  MATRIX  TO  BE  TRIANGULAR IZED 

C 

C  OUTPUT. . 

C 

C  A  CONTAINS  AN  UPPER  TRIANGULAR  MATRIX  U  AND  A  PREMUTED 
C  VERSION  OF  A  LOWER  TRIANGULAR  MATRIX  I~L  SO  THAT 

C  (PERMUTATION  MATRIX) *A=L*U 

C 

C  COND  =  AN  ESTIMATE  OF  THE  CONDITION  OF  A. 

C  FOR  THE  LINEAR  SYSTEM  A*X  -  B  ,  CHANGES  IN  A  AND  B 

C  MAY  CAUSE  CHANGES  COND  TIMES  AS  LARGE  IN  X. 

C  IF  COND+l.O  .EQ.  COND  ,  A  IS  SINGULAR  TO  WORKING 

C  PRECISION.  COND  IS  SET  TO  1.0D+32  IF  EXACT 

C  SINGULARITY  IS  DETECTED. 

C 

C  IPVT  -  THE  PIVOT  VECTOR 

C  IPVT (K)  -  THE  INDEX  OF  THE  K-TH  PIVOT  ROW 

C  IPVT (N)  -  (-1)** (NUMBER  OF  INTERCHANGES) 

C 

C  WORK  SPACE..  THE  VECTOR  WORK  MUST  BE  DECLARED  AND  INCLUDED 
C  IN  THE  CALL.  ITS  INPUT  CONTENTS  ARE  IGNORED. 

C  ITS  OUTPUT  CONTENTS  ARE  USUALLY  UNIMPORTANT. 


144 


no  n  noon 


THE  DETERMINANT  OF  A  CAN  BE  OBTAINED  ON  OUTPUT  BY 
DET(A)  -  IPVT(N)  *  A(1 , 1)  *  A(2,2)  *  ...  *  A(N,N)  . 


C 


C 

C 

C 


C 


C 

C 

C 

C 

C 

C 

C 

C 

C 


IPVT (N) =1 

IF  (N.EQ.l)  GO  TO  150 
NM1-N-1 

COMPUTE  THE  1-NORM  OF  A  . 

ANORM-O.DO 
DO  20  J-l.N 
T=0.D0 
DO  10  I«1,N 
10  T“T+DABS(A(I, J)) 

IF  (T.GT.ANORM)  ANORM=T 
20  CONTINUE 

DO  GAUSSIAN  ELIMINATION  WITH  PARTIAL 
PIVOTING. 

DO  70  K=1 ,NM1 
KPl*K+l 

FIND  THE  PIVOT. 

M=K 

DO  30  I=KP1 ,N 

IF  (DABS(A(I,K)) .GT.DABS(A(M,K)))  M=I 
30  CONTINUE 
IPVT (K) «M 

IF  (M.NE.K)  IPVT(N)=-IPVT(N) 

T-A(M,K) 

A(M,K)=A(K,K) 

A(K,K)=T 

SKIP  THE  ELIMINATION  STEP  IF  PIVOT  IS  ZERO. 
IF  (T.EQ.O.DO)  GO  TO  70 

COMPUTE  THE  MULTIPLIERS. 

DO  40  I=KP1,N 
40  A(I,K)=-A(I,K)/T 

INTERCHANGE  AND  ELIMINATE  BY  COLUMNS. 

DO  60  J=KP1,N 
T=A(M, J) 

A(M,  J)=A(K,  J) 

A(K,J)-T 

IF  (T.EQ.O.DO)  GO  TO  60 
DO  50  I»KP1,N 

50  A(I,  J)=A(I,  J)+A(I,K)*T 
60  CONTINUE 
70  CONTINUE 

COND  «  (1-NORM  OF  A)* (AN  ESTIMATE  OF  THE  1-NORM  OF  A- INVERSE) 

THE  ESTIMATE  IS  OBTAINED  BY  ONE  STEP  OF  INVERSE  ITERATION  FOR  THE 
SMALL  SINGULAR  VECTOR.  THIS  INVOLVES  SOLVING  TWO  SYSTEMS 
OF  EQUATIONS,  (A-TRANSPOSE) *Y  «  E  AND  A*Z  -  Y  WHERE  E 
IS  A  VECTOR  OF  +1  OR  -1  COMPONENTS  CHOSEN  TO  CAUSS  GROWTH  IN  Y. 
ESTIMATE  -  (1-NORM  OF  Z) / (1-NORM  OF  Y) 

SOLVE  (A-TRANSPOSE) *Y  -  E  . 
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DO  100  K**l  ,N 
T=0.D0 

IF  (K.EQ.l)  GO  TO  90 
Khl=K-l 
DO  80  I-l.KMl 
80  T=T+A(I,K)*WORK(I) 

90  EK-1.D0 

IF  (T.LT.O.DO)  EK=- 1 . DO 
IF  (A(K,K) .EQ.O.DO)  GO  TO  160 
100  WORK(K) =- (EK+T) /a(k,k) 

DO  120  KB=1 ,NM1 
K=N-KB 
T=O.DO 
KP1-K+1 

DO  110  I=KPl,N 
110  T=T+A(l,K)*WORK(K) 

WORK (K) =T 
M=IPVT(K) 

IF  (M.EQ.K)  GO  TO  120 
T=WORK (M) 

WORK (M) =WORK (K) 

WORK (K) =T 
120  CONTINUE 
C 

YNORM=0 . DO 
DO  130  1=1, N 

130  YNORM=YNORM+DABS(WORK(l)) 

C 

C  SOLVE  A*Z  =  Y 

CALL  SOLVE  (NDIM,N,A,WORK, IPVT) 

C 

ZNORM=0 . DO 
DO  140  1=1, N 

140  ZNORM=ZNORM+DABS (WORK (i)  ) 

C 

C  ESTIMATE  THE  CONDITION. 

COND=ANORM*ZNORM/ ynorm 
IF  (COND.LT. 1. DO)  C0ND=1 . DO 
RETURN 

C  1-BY-l  CASE.. 

150  C0ND= 1 . DO 

IF  (A ( 1 , 1 )  .NE.O.DO)  R£:'JRN 

EXACT  SINGULARITY 

160  C0ND=1 . 0D32 
RETURN 
END 

SUBROUTINE  SOLVE  (NDIM,N,A,B, IPVT) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(NDIM,N) , B (N) , IPVT (N) 

C 

C  SOLVES  A  LINEAR  SYSTEM,  A*X  =  B 

C  DO  NOT  SOLVE  THE  SYSTEM  IF  DECOMP  HAS  DETECTED  SINGULARITY. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


-COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS-,  BY  G.  E.  FORSYTHE, 
M.  A.  MALCOLM,  AND  C.  B.  MOLER  (PRENTICE-HALL,  1977) 

INPUT. . 

NDIM  -  DECLARED  ROW  DIMENSION  OF  ARRAY  CONTAINING  A 
N  -  ORDER  OF  MATRIX 

A  *=  TRIANGULARIZED  MATRIX  OBTAINED  FROM  SUBROUTINE  DECOMP 

B  =  RIGHT  HAND  SIDE  VECTOR 

IPVT  =  PIVOT  VECTOR  OBTAINED  FROM  DECOMP 

OUTPUT. . 

B  =  SOLUTION  VECTOR,  X 


DO  THE  FORWARD  ELIMINATION. 

IF  (N.EQ.l)  GO  TO  50 
NM1=N-1 
DO  20  K=1 ,NM1 
KP1=K+1 
M=IPVT(K) 

T=B  (M) 

B  (M)  =B  (K) 

B (K) -T 

DO  10  I=KP1,N 
10  B (I) =B (I) +A(I ,K) *T 
20  CONTINUE 

NOW  DO  THE  BACK  SUBSTITUTION. 

DO  AO  KB=1 ,NMl 
KM1=N-KB 
K=KM1+1 

B(K)=B  (K)  /A(K,K) 

T=-B (K) 

DO  30  1=1, KM1 
30  B(I)=B(I)+A(I,K)*T 

40  CONTINUE 
50  B(1)=B(1)/A(1,1) 

RETURN 

END 
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